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ABSTRACT 


A  general  purpose  software  package  was  developed  to 
perform  nonlinear  constrained  optimization  of  user-defined 
engineering  design  problems  of  significant  complexity  using 
desktop  computers.  The  package,  designated  Microcomputer- 
based  Design  Optimization  Tool  (MDOT),  will  accept  nonlinear 
functions  of  up  to  ten  variables,  which  may  be  bounded,  with 
as  many  as  fifty  constraints.  It  was  implemented  on  a 
Hewlett-Packard  Model  85  microcomputer  with  32  Kbytes  of 
random  access  memory. 

MDOT  employs  the  method  of  feasible  directions  for  con¬ 
strained  optimization,  and  a  variable  metric  method  for 

unconstrained  functions.  It  is  interactive,  provides  for 

* 

monitoring  the  optimization  progress,  and  can  be  interrupted 
to  restart  from  a  new  point  in  the  design  space.  Typical 
applications  of  MDOT  are  in  the  design  of  machine  com¬ 
ponents,  composite  laminates,  and  piping  systems. 
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NOMENCLATURE 


The  nomenclature  defined  here  is  that  used  in  the  text. 
Definitions  of  parameters  associated  with  the  program  code 
accompany  the  code  listings  in  Appendix  B.  Boldface  charac¬ 
ters  denote  vectors  or  matrices. 

a  One-dimensional  search  step  length 

D  Inverse  Hessian  approximation  update  matrix 

F  Objective  function 

G  Vector  of  inequality  constraints 

H  Hessian  matrix,  or  an  approximation  to  its  inverse 
1  Vector  of  design  variable  lower  bounds 
m  Number  of  inequality  constraints 

n  Number  of  design  variables 

p  Vector  used  in  constructing  D 

q  Iteration  number 

S  Search  direction  vector 

s  Scaler  used  in  constructing  D 

t  Scaler  used  in  constructing  D 

u  Vector  of  design  variable  upper  bounds 

w  Scaler  used  in  constructing  D 

X  Vector  of  design  variables 

y  Vector  used  in  constructing  D 

VF  objective  function  gradient  vector 
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I.  INTRODUCTION 


A.  OBJECTIVE 

This  thesis  presents,  and  describes  the  development  of, 
a  computer  software  package:  "Microcomputer-based  Design 
Optimization  Tool"  (MDOT).  The  motivation  for  this  work 
stemmed  from  the  lack  of  available  general  purpose  programs 
capable  of  performing  nonlinear  constrained  optimization  of 
engineering  design  problems  of  significant  complexity  using 
desktop  computers. 

In  a  more  general  3ense,  MDOT  is  intended  to  help  focus 
attention  on  the  versatility  and  computational  power  of 
microcomputers.  These  machines  are  a  potentially  very 
valuable  resource  which  is  just  beginning  to  be  tapped  by 
the  engineering  design  community. 

The  remainder  of  Chapter  I  is  devoted  to  an  overview  of 
where  microcomputers  stand  in  engineering  design  ,  where 
MDOT  stands  amid  the  optimization  software  currently  avail¬ 
able,  and  the  implementation  of  MDOT.  In  Chapter  II,  a 
general  description  of  optimization  concepts  and  methods, 
and  their  application  to  engineering  design,  is  presented. 
In  Chapter  III,  the  program  development  of  MDOT  is  de¬ 
scribed,  and  flowcharts  of  the  algorithms  coded  are  pro¬ 
vided.  In  Chapter  IV,  the  test  problems  which  were  used  to 
validate  MDOT  are  described,  along  with  the  solutions 
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obtained.  Chapter  V  is  a  brief  summary.  Appendix  A  is  the 
MDOT  user  manual.  Appendix  B  contains  an  annotated  listing 
of  the  MDOT  program  code. 

B.  OVERVIEW 

There  are  desktop  computers  available  today  with  memory 
size  and  computational  speed  in  excess  of  those  of  the 
mainframes  of  just  a  few  years  ago.  The  fact  that  their 
capabilities  are  not  yet  being  fully  exploited  in  the  day- 
to-day  process  of  engineering  design  can  be  attributed  in 
part  to  the  lack  of  available  software.  As  Falk  [Ref.  1: 
p.50]  observes,  "...engineers. ..have  little  time  or  pa¬ 
tience  to  do  computer  programming."  Even  among  those  engin¬ 
eers  who  have  the  time  and  patience,  there  persists  a  reluc¬ 
tance  to  program  on  microcomputers  because  of  a  perceived 
lack  of  general  purpose  utility  or  under-estimation  of  the 
capability  of  these  machines. 

Design  optimization  is  a  concept  which,  similar  to  the 
desktop  computer,  has  received  "mixed  reviews"  from  the 
engineering  disciplines.  While  there  are  few  who  would  ques¬ 
tion  the  virtue  of  seeking  the  "best"  solution  to  a  problem, 
there  are  many  who  are  reluctant  to  relinquish  to  a  computer 
what  they  see  as  the  engineers'  proprietary  decisions  in  the 
design  process  .  Our  whirlwind  courtship  of  computer  aided 
design  (CAD)  is  being  tempered  somewhat  by  a  counter  trend 
back  toward  "human  aided  design". 
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It  was  within  the  framework  of  these  two  ideas  that  the 
development  of  MDOT  was  undertaken;  not  only  to  make  avail¬ 
able  a  useful  interactive  design  optimization  program,  but 
to  demonstrate  that  a  powerful  general  purpose  problem  sol¬ 
ver  can  be  implemented  in  a  microcomputer.  Little  knowledge 
of  programming  is  required  of  the  MDOT  user.  Problem  entry 
and  program  execution  are  convenient.  The  interactive  fea¬ 
tures  of  the  code  permit  the  design  engineer  to  keep  in 
close  touch  with  the  progress  of  the  problem  solution  and  to 
interrupt  program  execution  to  make  parameter  adjustments 
based  on  engineering  judgement. 

The  chief  advantage  that  microcomputers  enjoy  over  main¬ 
frames  is  their  low  cost.  Small  computers  are  typically 
purchased  outright,  so  that  their  use  incurs  no  additional 
expense  for  connection  or  run  time.  Their  major  disadvantage 
is  comparatively  slow  computational  speed,  but  it  is  doubt¬ 
ful  that  engineering  design  ever  progresses  so  rapidly  as  to 
make  it  imperative  that  a  solution  be  obtained  in  seconds 
rather  than  minutes.  In  any  case,  this  disparity  is  rapidly 
disappearing  with  the  development  of  16  and  32-bit  micro¬ 
processor-based  desktop  computers  [Ref.  2:  p.2]. 

C.  OPTIMIZATION  SOFTWARE  CURRENTLY  AVAILABLE 

There  are  several  powerful  general  purpose  optimization 
programs  available,  such  as  COPES/CONMIN  [Ref.  3],  which  can 
deal  with  a  wide  range  of  design  problems.  These  programs 
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must  reside  in  a  mainframe  computer,  and  their  use  can  be 
cumbersome,  especially  for  the  occasional  user.  At  the  other 
extreme  are  those  codes  developed  for  use  in  computers  with 
limited  memory.  Typically  these  are  special  purpose  programs 
employing  zero  order  or  simple  first  order  methods,  such  as 
random  search  or  steepest  descent,  capable  of  handling  only 
relatively  small  problems.  They  are  convenient,  but  of  lim¬ 
ited  usefulness. 

The  gap  between  these  two  categories  requires  that  opti¬ 
mization  of  the  great  number  of  general  design  and  analysis 
problems  which  are  on  a  scale  that  could  easily  be  handled 
by  small  computers  be  done  on  a  mainframe  or  not  at  all. 
This  often  leads  to  overmodeling,  wherein  a  relatively 
simple  problem  is  unnecessarily  made  more  complicated  in 
order  to  more  fully  utilize  the  machine  capability  or  to 
justify  the  expense  of  computer  services.  MDOT  was  developed 
specifically  to  bridge  this  gap. 

MDOT  provides  the  design  engineer  with  a  convenient  tool 
for  optimization  of  nonlinear  problems  in  up  to  ten  bounded 
independent  variables  subject  to  as  many  as  fifty  inequality 
constraints.  All  that  is  required  is  access  to  a  desktop 
computer,  and  today  there  are  certainly  few  engineers  who 


lack  this 


D.  IMPLEMENTATION  OF  MDOT 


The  program  development  for  MDOT  was  done  on  a  Hewlett- 
Packard  model  85A  microcomputer,  which  is  built  around  an  8- 
bit  microprocessor.  This  particular  machine  is  so  often  used 
as  a  data  aquisition  system  controller  that  its  stand-alone 
computational  capability  may  frequently  be  overlooked.  A 
versatile,  engineering  oriented  computer  with  high  machine 
precision,  it  is  nonetheless  on  the  low  end  of  the  memory 
size  scale  with  just  16  Kbytes.  The  computer’s  capabilities 
were  enhanced  by  the  addition  of  a  16  Kbyte  memory  extension 
module  and  three  read-only-memory  modules:  the  matrix,  ad¬ 
vanced  programming,  and  printer/plotter  ROMs.  As  configured, 
there  were  just  over  30  Kbytes  of  memory  available  for 
programming. 

MDOT  was  written  in  HPBASIC,  which  differs  in  some  re¬ 
spects  from  standard  BASIC.  As  such,  in  its  present  form 
MDOT  is  limited  to  use  in  the  HP  series  80  computers.  With¬ 
out  much  difficulty,  though,  the  code  could  be  translated 
and  run  on  almost  any  available  hardware.  Additional  com¬ 
ments  concerning  transferrability  of  the  code  are  presented 
in  chapter  III. 


OPTIMIZATION 


A.  APPLICATION  TO  ENGINEERING  DESIGN 


Design  has  been  described  as  the  creative  process 
through  which  the  engineering  profession  develops  devices, 
processes  and  systems  to  fill  the  needs  of  man  [Ref.  M: 
p.170].  A  "need”  must  first  be  defined  in  terms  of  specific 
requirements  which  the  design  must  meet.  Then  the  engineer, 
drawing  on  available  resources,  synthesizes  proposed  solu¬ 
tions  to  meet  these  requirements.  Many  such  designs  may,  and 
usually  do,  exist.  Thus  there  arises  the  subproblem  of 
finding  the  best  of  these  designs,  and  the  inherently  iter- 
ative  nature  of  design. 

Traditionally,  the  solution  to  this  subproblem  was 
sought  through  comparitive  analysis  of  a  reasonable  number 
of  alternative  designs;  a  tedious  and  expensive  procedure 
for  problems  of  even  moderate  complexity.  The  recent  devel¬ 
opment  of  a  broad  range  of  very  useful  CAD  software  has  made 
it  possible  to  remove  a  good  deal  of  the  tedium,  and,  per¬ 
haps  to  a  lesser  degree,  the  expense  of  engineering  design. 
For  the  most  part,  though,  these  tools  have  made  no  fun¬ 
damental  change  in  the  approach  taken  to  solve  the  design 
problem.  What  they  have  done  is  redefined  the  phrase  "a 
reasonable  number  of  alternative  designs".  By  programming 


the  analysis  and  comparison  tasks  into  a  computer,  the 
engineer  is  able  to  consider  many  more  possible  solutions  in 
the  same  amount  of  time. 

Optimization  methods  are  a  significant  extension  of  the 
CAD  concept  in  that  they  enable  the  engineer  to  exploit  the 
capabilities  of  the  computer  over  the  entire  scope  of  the 
design  process.  In  optimization,  the  computer  is  tasked  not 
only  with  analysis  and  comparison  of  previously  selected 
designs,  but  with  selection  of  the  designs  to  be  considered 
in  subsequent  iterations  as  well.  Since  this  intermediate 
design  selection  can  be  quite  complex,  closing  the  design 
loop  in  the  computer  can  lead  to  a  considerable  savings  of 
time  and  effort  in  the  search  for  the  optimum. 

B.  THE  NATURE  OF  THE  PROBLEM 

Fundamental  to  the  economical  solution  of  the  design 
prohlem  is  that  it  be  quantified  and  formulated  mathematic¬ 
ally  to  permit  conceptual,  rather  than  physical,  manipul¬ 
ation  of  resources.  The  design,  then,  is  specified  by  ass¬ 
igning  values  to  a  set  of  independent  variables  which  repre¬ 
sent  its  physical  characteristics.  The  measure  of  goodness 
of  the  design  used  for  comparison  with  others  is  expressed 
as  some  functional  relationship  between  these  variables.  The 
requirements  placed  on  the  design,  as  well  as  the  physical 
limitations  of  the  design  itself,  define  a  region  in  the 
multi-dimensional  mathematical  design  space.  The  design  must 
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reside  inside  this  region  to  be  acceptable.  In  the  termin¬ 
ology  of  optimization,  the  measure  of  goodness  is  the 
objective  function,  the  requirements  and  limitations  are 
constraints,  and  designs  which  fall  within  the  region 
bounded  by  the  constraints  are  considered  feasible.  As 
examples,  the  set  of  independent  design  variables  might  be 
the  cross-sectional  dimensions  of  a  structural  element,  the 
objective  function  its  weight,  and  the  constraints  its 
maximum  allowable  stress  and  size  limitations. 

In  general,  the  formulation  of  a  design  problem  leads  to 
an  objective  and  a  number  of  constraints,  all  of  which  may 
be  linear  or  nonlinear  functions,  explicit  or  implicit  in 
many  design  variables  which  themselves  are  subject  to  limi¬ 
tations,  called  bounds  or  side  constraints.  Stated  mathe¬ 
matically  [Ref.  5:  p-9],  the  design  optimization  problem  is 
to 


Minimize: 

F(X) 

objective  function 

Subject  to: 

Gj ( X) <0 , 

j  =  1,m 

inequality  constraints 

1 ^ 

i  =  1  ,n 

side  constraints 

Where : 

XT={X1,X2, 

c 

X 

• 

• 

• 

design  variables 

1T={11,12,. 

lower  bounds  on  X 

uT={u1 ,u2, 

•  •  •  Ujj  } 

upper  bounds  on  X 

If  in  the  problem  formulation  it  is  more  convenient  to 
define  the  objective  function  as  a  quantity  which  is  to  be 
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maximized,  such  as  efficiency  or  torque,  then  the  above 
statement  may  simply  be  modified  to  read  "Minimize:  -F(X)". 

C.  THE  NATURE  OF  THE  SOLUTION 

Optimization  is  an  application  of  mathematical  theory 
concerning  identification  of  the  extrema  of  functions.  In 
multivariable  calculus,  for  example,  the  method  of  Lagrange 
multipliers  is  developed,  which  provides  a  closed  form 
solution  for  the  extremum  of  a  constrained  function.  While 
useful  for  demonstrating  concepts  and  developing  methods, 
such  analytical  techniques  are  not  practical  for  solving  any 
but  the  simplest  of  problems.  Design  optimization  methods 
involve  numerical  approximation  techniques  and  iterative 
search  schemes.  They  are  ideally  suited  to,  and  in  fact  made 
practical  only  through,  the  use  of  digital  computers. 

Many  optimization  algorithms  have  been  developed  around 
widely  varying  strategies.  Common  to  most  are  the  three 
basic  tasks  that  make  up  one  iteration  of  the  solution  loop: 

1.  Selection  of  a  direction  in  the  design  space  along 

which  to  search. 

2.  A  search  for  the  most  improved  design  in  this 
direction. 

3.  Convergence  testing  to  determine  when  the  optimum 
design  has  been  found. 


For  unconstrained  problems,  these  tasks  are  relatively 
straightforward.  Addition  of  a  constraint  set  may,  depending 
on  the  sophistication  of  the  method  employed,  complicate  the 
first  two  steps  considerably. 

Except  in  the  case  of  zero  order  methods,  selection  of  a 
search  direction  involves  calculation  of  partial  deriv¬ 
atives,  for  which  general  purpose  optimizers  use  numerical 
techniques,  such  as  finite  forward  differences.  At  any  point 
in  the  design  space,  the  negative  of  the  gradient  of  the 
objective  function  indicates  the  direction  in  which  the 
objective  function  is  most  rapidly  decreasing.  This  may  not 
be  the  best  direction  in  which  to  search,  however,  if  the 
objective  function  is  highly  nonlinear  or  if  the  design  is 
near  one  or  more  constraints.  Efficient  algorithms  variously 
employ  constraint  gradients,  Hessian  matrix  approximations, 
and  previous  iteration  history  information  in  addition  to 
objective  function  gradients  to  select  the  search  direction. 

Finding  the  best  improved  design  along  a  line  in  the 
specified  direction  is  termed  a  "one-dimensional  search", 
because  the  objective  and  constraints  are  treated  as  func¬ 
tions  only  of  the  "distance"  along  this  line  from  the  cur¬ 
rent  design  point.  Techniques  employed  in  the  one-dimen¬ 
sional  search  include  the  golden  section  and  Fibonacci 
methods,  polynomial  approximations,  and  combinations  of 
these  [Refs.  5,6],  If  constraints  are  present,  the  best 


violated,  the  design  is  infeasible,  so  the  search  algorithm 
must  seek  the  point  at  which  the  objective  function  is 
minimized  while  remaining  inside  the  feasible  region. 

Part  of  the  optimization  problem  formulation  is  specifi¬ 
cation  of  an  initial  design  point  from  which  to  start  the 
solution  process.  For  constrained  problems,  the  possibility 
exists  that  this  initial  design  will  be  infeasible.  To 
provide  for  this,  the  search  direction  routine  must  find  the 
direction  which  will  yield  the  shortest  path  to  the  feasible 
region,  and  the  one-dimensional  search  algorithm  must  allow 
for  the  possible  necessity  of  increasing  the  objective  func¬ 
tion  in  order  to  attain  feasibility. 

Convergence  to  a  global  optimum  generally  cannot  be 
guaranteed.  Theory  provides  the  Kuhn-Tucker  conditions 
necessary  for  the  existence  of  an  optimum,  but  these  are 
neither  convenient  to  evaluate  nor  sufficient  to  define 
optimality  [Ref.  5:  pp.  17-20],  In  practice,  convergence  is 
typically  considered  to  be  indicated  by  one  or  more  of  the 
following: 

1.  Failure  to  find  a  search  direction  which  will  lead  to 
an  improved  design. 

2.  Given  a  direction,  failure  to  find  any  significant 
search  step  length  to  improve  the  design. 

3.  Finding  no  appreciable  design  improvement  over  a 
specified  number  of  iterations. 


If  the  possibility  of  local  minima  exists,  the  optimi¬ 
zation  should  be  repeated  from  several  different  initial 
designs.  For  constrained  problems,  the  optimization  process 
may  fail  to  find  any  feasible  solution,  in  which  case  the 
problem  must  be  reformulated. 

An  iteration  in  the  optimization  solution,  then,  may  be 
summarized.  Beginning  from  the  current  design  point  Xq,  a 
search  direction,  S3,  is  determined.  Then  the  one-dimen¬ 
sional  search  is  conducted  to  find  the  "distance",  a,  along 
Sq,  which  yields  the  best  improved  design.  The  design  is 
then  updated  as 

X3+1  =  X3  +  as3 

at  which  point  the  objective  function  is  reevaluated  and  the 
design  checked  for  convergence. 

Optimization  algorithmic  efficiency  and  convergence 
behavior  are  affected  by  the  mathematical  characteristics  of 
the  problem.  As  numerical  methods,  they  are  susceptible  to 
ill-conditioning.  Truncation  and  round-off  errors,  which  are 
an  unavoidable  consequence  of  the  use  of  digital  computers, 
aggravate  this.  Given  an  optimizer  suitable  to  the  problem 
type,  careful  problem  formulation,  as  discussed  in  Appendix 
A,  is  the  best  insurance  against  poor  optimizer  performance. 


III.  PROGRAM  DEVELOPMENT 


A.  BASIC  CRITERIA 

At  the  outset  of  the  program  development,  four  basic 
criteria  were  established  to  be  met  by  MDOT:  utility,  mini¬ 
mization  of  required  memory,  user  convenience,  and  reduction 
of  problem  run  time.  At  points  where  conflict  existed 
between  them,  these  criteria  were  prioritized  in  the  order 
listed.  Few  such  compromises  were  necessary,  as  the  require¬ 
ments  were  found  to  be  generally  complimentary. 

The  utility  criterion  meant  that  MDOT  should  be  a 
general  purpose  optimizer  which  could  be  applied  to  a  wide 
range  of  design  problems.  Minimization  of  memory  was 
dictated  by  the  limitaions  of  microcomputers,  and  affected 
not  only  algorithm  selection  and  problem  size,  but  many 
aspects  of  the  actual  coding  as  well.  User  convenience 
considerations  drove  the  development  of  those  portions  of 
the  code  which  are  interactive,  and  those  involved  with 
problem  entry  and  output  options.  Reduction  of  problem  run 
time  was  a  factor  throughout  the  development,  most  notably 
in  the  incorporation  of  optimization  progress  display  and 
the  interrupt/restart  option. 


B.  MDOT  ALGORITHMS 


1.  Algorithm  Selection 

Of  the  many  optimization  algorithms  available,  the 
zero  order  methods,  as  well  as  the  simpler  of  the  first 
order  methods,  were  ruled  out  on  the  basis  of  their  lack  of 
general  purpose  utility.  Others,  including  linear  and  quad¬ 
ratic  programming  types,  were  eliminated  because  of  their 
excessive  memory  requirements  [Refs.  7,83.  Finally,  the  need 
to  reduce  problem  run  time  while  retaining  utility  lead  to 
the  selection  of  two  algorithms,  each  capable  of  nonlinear 
multidimensional  optimization.  The  first  is  a  variable 
metric  method  which  is  used  in  MDOT  for  unconstrained  opti¬ 
mization,  the  second  is  a  method  of  feasible  directions,  for 
minimizing  constrained  functions. 

The  selection  of  the  one-dimensional  search  strategy 
to  be  employed  was  driven  by  the  need  to  reduce  problem  run 
time.  There  is  a  trade-off  to  be  made  between  the  precision 
to  which  the  search  step  length  determination  is  made  and 
the  time  required  for  each  optimization  iteration.  Both  the 
golden  section  and  Fibonacci  search  methods  can  attain  very 
precise  solutions,  but  to  do  so  they  become  computationally 
expensive.  In  MDOT  the  one-dimensional  search  routines  were 
designed  to  seek  a  less  precise  step  length  solution  in 
order  to  complete  each  iteration  more  quickly.  The  method 
employed  in  both  the  optimizers  estimates  an  initial  step 
length  based  on  a  reasonable  change  in  objective  function 
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magnitude.  The  golden  section  ratio  is  then  used  to  estab¬ 
lish  bounds  on  the  solution,  which  is  finally  refined  by 
polynomial  approximation. 

2.  The  Unconstrained  Optimizer 

From  the  derivation  of  the  Kuhn-Tucker  conaitions, 
it  is  known  that  if  at  some  point  X*  the  objective  function 
F(X*)  has  a  local  minimum,  then  the  gradient  of  the  object¬ 
ive  at  this  point,  7F(X*),  must  vanish  and  the  Hessian 
matrix  H  must  be  positive  definite.  Combined  with  a  second 
order  Taylor  series  expansion  of  F(X)  about  some  point,  say 
X°,  near  the  minimum,  these  conditions  lead  to  an  expression 
for  the  direction  from  X°  to  X*  as 

X*  -  X°  s  -H(X°)"Wf(X°) 

In  practice,  determination  of  the  Hessian  matrix  by  finite 
difference  approximation,  as  well  as  inversion  of  the 
matrix,  would  be  computationally  so  expensive  as  to  outweigh 
the  theoretical  gain  in  algorithmic  efficiency. 

In  variable  metric  methods,  information  gathered  as 
the  optimization  progresses  is  used  to  develop  an  approxi¬ 
mation  to  the  inverse  of  the  Hessian  matrix,  which  is  then 
used  in  determining  the  search  direction.  As  such,  these 
first  order  methods  have  some  convergence  character ist ics 
comparable  to  those  of  second  order  methods.  The  algorithm 
for  the  variable  metric  method  is  shown  in  Fig.  1. 
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The  inverse  Hessian  approximation  is  initialized  as 
an  n  x  n  identity  matrix.  To  begin  each  iteration,  the 


search  direction  is  defined  as 

S9  =  -HVFCX^) 

After  the  one-dimensional  search  and  design  update,  H  is 

modified  as  , 

H^+  = 

where  the  form  of  the  update  matrix  D  determines  which  one 
of  a  family  of  variable  metric  methods  is  being  used.  In 
general,  D  is  defined  as 

Dq  =  3-+^J'  ppT  +  Hqy(Hqy)T  -  -[H^ypT+p(Hqy)T] 

s2  t  s 

where:  p  -  -  X^”1 

y  s  VF(X^)  -'7F(X^-1) 

s  =  p-y 

t  =  yTH^y 


Two  forms  of  D,  and  thus  two  variable  metric 
methods,  are  available  in  MDOT.  The  first  is  the  Davidon- 
Fletcher-Powell  method,  where  w  is  set  equal  to  zero.  The 
second  is  the  Broydon-Fletcher-Goldfarb-Shanno  method,  with 
w  equal  to  one  [Ref.  5:  pp.  92,93]-  As  the  convergence 
behavior  of  a  given  algorithm  can  be  somewhat  problem  depen¬ 
dent,  this  feature  allows  the  MDOT  user  to  compare  the 
results  of  two  variations  of  unconstrained  optimization. 
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The  one-dimensional  search  routine  employed  by  MDOT 
in  the  variable  metric  optimizer  first  finds  bounds  on  the 
unconstrained  minimum  of  the  objective  function,  then 
refines  the  minimum  by  a  three-point  cubic  polynomial 
approximation. 

3.  The  Constrained  Optimizer 

The  addition  of  a  set  of  constraints  to  the  optimiz¬ 
ation  problem  requires  that  more  sophisticated  techniques  be 
applied  to  its  solution,  particularly  in  the  determination 
of  a  search  direction  and  the  subsequent  one-dimensional 
search.  As  is  the  case  for  unconstrained  optimizers,  it  is 
generally  the  method  used  to  find  a  search  direction  which 
distinguishes  the  different  constrained  optimization  algo¬ 
rithms.  In  the  method  of  feasible  directions,  a  search 
direction  in  which  a  finite  step  will  reduce  the  objective 
function  is  termed  useable,  while  one  which  will  avoid 
constraint  violation  is  called  feasible.  The  direction  find¬ 
ing  problem  is  then  formulated  as  a  sub-optimization  task  to 
determine  the  best  of  the  possible  useable-feasible  direct¬ 
ions.'  MDOT  employs  the  algorithm  presented  by  Vanderplaats 
[Ref.  5:  pp.  163-1703  for  the  solution  of  this  sub-problem. 

As  shown  in  the  flowchart  of  Fig.  2,  the  feasible 
directions  optimizer  begins  as  a  simple  steepest  descent 
algorithm,  provided  the  initial  design  is  feasible.  Opti¬ 
mization  thus  proceeds  quickly  to  a  point  where  one  or  more 


constraints  are  encountered 


E: 


On  subsequent  iterations  the  sub-optimization  rou¬ 
tine  is  used  to  determine  search  directions  so  as  to  satisfy 
the  fundamental  requirement  that  the  optimum  design  be 
feasible.  The  constrained  optimizer  in  MDOT  will  accept 
initially  infeasible  designs.  In  this  case,  the  direction 
and  search  routines  are  modified  so  as  to  attain  feasibility 
as  quickly  as  possible.  Thus  a  feasible  direction  and  step 
length  are  sought  which  will  overcome  the  constraint  vio¬ 
lations,  even  at  the  expense  of  increasing  the  objective 
function.  Once  inside  the  feasible  region,  optimization 
proceeds  as  before. 

As  in  the  unconstrained  case,  in  the  one-dimensional 
search  routine  employed  in  MDOT  for  constrained  optimiz¬ 
ation,  bounds  on  the  solution  are  first  established,  fol¬ 
lowed  by  refinement  by  polynomial  approximation.  Here,  how¬ 
ever,  the  search  must  be  conducted  for  the  zeros  of  the 
constraint  functions  as  well  as  for  the  minimum  of  the 
objective.  The  step  length  selected  is  then  the  one  which 
yields  the  best  feasible  design.  Provision  must  also  be  made 
to  ensure  the  design  variables  remain  within  their  bounds 
(side  constraints).  In  MDOT,  if  at  any  time  during  the  one¬ 
dimensional  search  a  design  variable  is  found  to  exceed  an 
upper  or  lower  bound,  it  is  set  equal  to  the  value  of  the 
violated  bound. 
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The  relationships  between  the  modules  of  MDOT  are 
depicted  in  Figs.  3  and  4  for  unconstrained  and  constrained 
optimization,  respectively.  The  main  program  is  named 
"Autost"  because  this  signals  the  HP-85  operating  system  to 
load  and  run  this  program  automatically  when  the  computer  is 
powered  up  with  the  mass  storage  cartridge  inserted  in  the 
tape  drive.  All  the  other  modules  of  MDOT  are  subprograms 
which  are  called  into  memory  and  executed  as  needed  by 
Autost  or  another  subprogram.  Once  entered  into  main  memory, 
a  subprogram  resides  there  for  the  duration  of  the  optimiz¬ 
ation  unless  a  "SCRATCHSUB"  instruction  is  executed. 

Following  is  a  brief  description  of  the  function  of 
each  of  MDOT's  program  segments: 

Autost  MDOT  main  calling  program 

LOGO  Displays  introductory  (welcome)  graphic 

DEFALT  Sets  program  parameters  to  default  values 

PROB  Problem  entry,  evaluation  of  F  and  G 

CCONT  Control  of  constrained  optimization 

UCONT  Control  of  unconstrained  optimization 

ACON  Identification  of  active/violated  constraints 

GRAD  Calculation  of  gradients  of  F  and  G 

DIRECT  Direction  finding  subproblem  solver 

FDSRCH  Constrained  one-dimensional  search 

VMSRCH  Unconstrained  one-dimensional  search 

Update  approximation  to 


NEWH 


CONV 


Convergence  testing 
PROG  Optimization  progress  information 

TERM  Output  of  results  of  the  optimization 

The  PROB  subprogram  is  created  by  the  user  by 
editing  a  skeleton  problem  entry  code  which  is  stored  on  the 
tape.  The  edited  version  is  then  renamed  and  stored.  Autost 
queries  the  user  for  the  problem  name,  which  is  then  common 
to  all  subprograms.  Autost  makes  the  first  call  to  PROB  in 
order  to  select  the  appropriate  optimization  control  rou¬ 
tine.  Thereafter,  PROB  is  called  any  time  an  objective 
function  or  constraint  evaluation  is  required.  Both  LOGO  and 
DEFALT  are  called  by  Autost.  LOGO  generates  a  simple 

"welcome"  graphics  display,  and  is  scratched  from  memory 

* 

upon  execution.  DEFALT  initializes  a  number  of  program  para¬ 
meters  to  their  default  values,  as  defined  in  Appendix  B. 

PROG  is  called  upon  completion  of  each  iteration. 
This  subprogram  generates  the  user  selected  optimization 
progress  indicators.  Options  include  data  and  graphics 
displays  and  printed  output.  Based  on  the  progress  inform¬ 
ation  provided,  the  user  may  elect  to  continue  the  optimiz¬ 
ation,  restart  MDOT  from  a  different  initial  design,  or  stop 
and  reformulate  the  problem.  CONV  is  also  called  upon 
completion  of  each  iteration,  to  determine,  based  on  the 
convergence  criteria  set  by  DEFALT,  whether  the  optimum 
design  has  been  found. 
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TERM  is  called  to  end  the  optimization  and  generates 
the  output  of  results.  Termination  may  be  invoked  by  a 
number  of  conditions  in  other  than  meeting  convergence 
criteria.  If  the  optimizer  exceeds  a  specified  number  of 
iterations,  if  the  components  of  the  search  direction  vector 
are  all  essentially  zero,  if  no  search  step  length  can  be 
found  to  improve  the  design,  or  if  there  are  an  excessive 
number  of  violated  constraints,  MDOT  will  terminate.  In  any 
case,  TERM  will  generate  an  output  message  to  indicate  the 
condition  upon  which  the  decision  to  terminate  was  based, 
and  offer  the  user  the  option  of  editing  and  restarting  the 
program. 

C.  ADAPTATION  OF  MDOT  TO  OTHER  SYSTEMS 

Since  MDOT  is  coded  in  HPBASIC,  it  is  not  immediately 
transf errable  to  hardware  other  than  the  Hewlett-Packard 
series  80  desktop  computers.  Translation  of  the  package, 
either  into  another  version  of  BASIC  or  into  FORTRAN,  is 
certainly  a  "do-able"  project  which  would  significantly 
expand  the  applicability  of  MDOT. 

This  section  highlights  those  features  of  HPBASIC  used 
in  MDOT  which  would  have  the  greatest  impact  on  this  pro¬ 
ject.  They  are:  variable  name  assignment,  SUBPROGRAMS, 
matrix  manipulations,  and  graphics. 

The  limitation  of  HPBASIC  which  most  decreases  the  read¬ 
ability  of  the  code  is  that  variable  names  are  restricted  to 
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either  one  letter  or  a  letter  followed  by  a  digit.  One 
result  of  this  is  that  arrays  of  subscripted  variables  are 
sometimes  used  where  individual  characteristic  names  might 
otherwise  be  assigned.  A  feature  of  the  language  which  helps 
considerably,  not  only  to  overcome  this  limitation,  but  in 
programming  complicated  algorithms,  is  the  SUBPROGRAM. 

Similar  to  a  SUBROUTINE  in  FORTRAN,  the  SUBPROGRAM  is 
called  when  needed  and  variables  may  be  passed  either  by 
name  or  by  value.  In  HPBASIC  this  allows  for  the  use  of  the 
same  variable  name  to  denote  different  parameters  in  sep¬ 
arate  program  segments.  MDOT  makes  extensive  use  of  SUBPRO¬ 
GRAMS.  This  feature  is  not  available  in  all  versions  of 
BASIC.  Without  it,  the  translation  of  MDOT  would  be  more 
difficult,  but  still  possible  through  the  use  of  functions 
and  subroutines,  particularly  if  multi-character  variable 
names  are  permitted. 

Matrix  manipulation  i3  convenient  in  HPBASIC.  Operations 
such  as  matrix  multiplications,  transpositions,  dot  pro¬ 
ducts,  and  identifying  extreme  array  elements  are  all 
accomplished  through  simple  "MAT"  statements.  This  feature 
is  available  in  some  of  the  other  versions  of  BASIC,  but  not 
in  FORTRAN.  Without  it,  additional  subprogramming  would  be 
required  to  perform  these  operations. 

Graphics  capabilities  vary  widely  from  one  hardware 
manufacturer  to  another,  as  do  the  coding  instructions  used 
to  execute  the  displays.  It  is  likely  that  the  graphics  pro- 


grammed  into  MOOT  would  require  major  modification  to  make 
them  transferable. 

Two  more  details  of  the  HPBASIC  code  are  worthy  of  note: 
program  flags  and  the  g  symbol.  Program  flags  are  built-in 
indicators  which  can  be  set  to  1  or  cleared  to  0,  and  are 
used  in  MDOT  for  conditional  branching  decisions.  They  could 
easily  be  replaced  by  integer  variables.  The  g  symbol  is 
used  to  condense  the  code  and  thus  conserve  memory.  It 
simply  separates  multiple  executable  statements  on  one  pro¬ 
gram  line.  Without  this  feature,  each  statement  must  have 
its  own  line  number. 

D.  POTENTIAL  FOR  FUTURE  GROWTH 

Besides  its  obvious  potential  for  expanded  problem  size 
if  implemented  in  a  computer  with  a  larger  memory,  there  are 
many  refinements  and  additions  which  could  be  incorporated 
into  MDOT,  either  to  enhance  its  general  purpose  utility  or 
to  tailor  it  to  a  particular  type  of  problem.  Some  modifi¬ 
cations  for  improved  utility  might  involve  coding  additional 
algorithms,  automatic  design  variable  scaling,  and  the  hand¬ 
ling  of  equality  constraints. 

MDOT  could  be  customized  by  modification  of  the  problem 
input  subprogram,  graphics  display,  output  format,  or  the 
algorithms  themselves.  Coupling  of  MDOT  to  an  external  CAD 
or  analysis  code  also  presents  many  interesting  and  poten¬ 
tially  useful  possibilities.  One  such  configuration  might  be 
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to  use  MDOT  as  a  subprogram  called  to  perform  optimization 
on  localized  aspects  of  a  large  scale  design  problem  in 
conjunction  with  a  desktop  computer  CAD  system  [Ref.  91. 


IV.  TEST  CASES 


Validation  of  an  optimization  program  typically  consists 
of  testing  it  on  a  battery  of  representative  problems  to 
which  the  solutions  are  known  [Ref.  10].  Based  on  the 
results  of  such  tests,  a  number  of  yardsticks  exist  by  which 
the  optimizer  is  judged  relative  to  others.  These  can  be 
grouped  into  three  categories:  stability,  robustness,  and 
efficiency  [Ref.  8:  p.75]. 

An  optimizer  is  stable  if,  once  a  feasible  design  is 
attained,  the  objective  function  remains  non-increasing 
until  the  optimum  has  been  found.  A  robust  optimizer  is  one 
which  yields  a  valid  solution  given  a  poor  initial  approxi¬ 
mation.  Efficiency  refers  either  to  the  number  of  function 
and  derivative  evaluations  required  in  the  solution  or  to 
the  problem  run  time.  These  two  measures  of  efficiency  are 
closely  related  if  comparing  different  optimizers  run  on  the 
same  machine,  since  function  and  derivative  evaluations  are 
typically  costly  operations. 

To  these  measures  of  an  optimizer’s  performance,  a 
fourth  category  should  be  added;  that  of  utility.  Given  a 
stable  and  robust  optimizer,  there  are  characteristics  in 
addition  to  its  efficiency  which  should  be  considered  in 
determining  its  utility.  Problem  size  and  type  solvable  by 
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the  program  are  of  fundamental  importance,  as  is  memory 
storage  required  to  run  it.  There  are  trade-offs  involved 
here  in  program  development  between  the  sophistication  of 
the  algorithms  used,  speed  of  convergence,  hardware 
capabilities,  cost  of  run  time,  and  user  convenience.  The 
utility  of  the  optimizer  is  an  indication  of  how  these 
trade-offs  were  made,  and  involves  much  more  than  just 
efficiency. 

A.  UNCONSTRAINED  TEST  PROBLEM 

Among  the  unconstrained  test  problems  run  on  MDOT  was 
the  so-called  "banana"  function: 

F(X)  s  lOX^  -  20X12X2  +  10X22  +  X^  -  2X^  +  5 

which  has  an  optimum  of  F(1,1)  =  4.0.  This  function  derives 
its  name  from  the  shape  of  the  contours  of  constant  objec¬ 
tive  function  (Fig.  5).  Although  only  two-dimensional,  the 
banana  function  is  a  good  test  of  an  unconstrained  optimizer 
because  the  objective  function  surface  becomes  a  steep, 
narrow,  curved  "valley"  as  the  optimum  is  approached.  An 
inefficient  optimizer  will  tend  to  "zig-zag"  in  such  a 
design  space,  resulting  in  slow  convergence  near  the  solu¬ 
tion,  while  a  non-robust  optimizer  will  tend  to  terminate 
prematurely. 
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Results  of  the  performance  of  MDOT  on  the  banana 


function  are  summarized  below. 

=  -1.0 
X2°  =  1.5 
F°  r  10.5 

X/  =  0.95979 
X2*  =  0.91442 
F*  =  4.002 

MDOT  arrived  at  this  solution  in  about  one  minute,  after  12 
iterations  and  66  function  evaluations. 

B.  CONSTRAINED  TEST  PROBLEM 

Among  the  constrained  test  problems  run  on  MDOT  was  the 
cantilevered  beam  problem  posed  by  Vanderplaats  [Ref.  3: 
p.8],  as  illustrated  in  Fig.  6.  The  objective  function  in 
this  case  is  the  volume  of  the  beam,  for  which  a  theoretical 
optimum  of  6603.9  is  known.  Results  of  the  performance  of 
MDOT  on  this  beam  design  problem  are  summarized  below. 


Initial  design: 


Optimum: 


Initial  design: 


X!°  =  3.5 
X2°  =  16.0 
F°  =  11200 


v  >VV1 


Cant i  1  evened  Beam: 


P-1Q.  □□□  lb 


1-200  in 


Design  Variables: 
Objective  Function: 
Constraints: 


Where: 


b,  h 
Volume 

Bending  Stress  <  20,000  psi 
Deflection  <  1.0  in 
Ratio  of  h  to  b  <  10.  0 
0.5  <  b  <  5.0 
1.0  <  h  <  20.  0 

Bending  Stress  =  Me/ I 
Deflection  =  P1V3EI 


Fig.  6  Constrained  Test  Problem 


rV 


Optimum: 


1.8261 


X2*  =  18.174 

F*  =  6637.5 

G/  =  -0.024 
G2*  =  -0.00656 
G3*  =  -0.00522 

MDOT  arrived  at  this  solution  in  about  two  minutes,  after  11 
iterations  and  56  function  evaluations. 

With  refinement  of  the  algorithms,  improvement  could 
likely  be  realized  in  the  program  performance.  All  variables 
in  MDOT  are  declared  "SHORT",  which  in  HPBASIC  means  they 
are  carried  to  5  digits."  In  a  machine  with  just  64  Kbytes  of 
memory,  this  could  be  changed  to  "REAL",  in  which  case  9 
digits  would  be  carried,  with  an  attendant  improvement  in 
the  precision  of  the  solution.  The  number  of  function  eval¬ 
uations,  and  thus  the  problem  run  time,  could  be  decreased 
by  modification  of  the  algorithms  such  that  gradients  are 
not  calculated  in  every  iteration.  Also,  the  efficiency  of 
the  one-dimensional  search  routines  could  be  improved  by 
distinguishing  between  linear  and  nonlinear  constraint 
functions. 
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V.  SUMMARY 


Optimization  is  a  useful  tool  in  engineering  design.  The 
desktop  computer  is  the  vehicle  through  which  this  tool  can 
be  made  widely  available,  convenient,  and  inexpensive.  The 
development  of  MDOT  affirms  the  feasibility  of  implementing 
a  powerful  general  purpose  optimization  algorithm  in  a  com¬ 
puter  with  limited  memory. 

The  applicability  of  MDOT  could  be  expanded  through  con¬ 
version  to  standard  BASIC  or  translation  to  FORTRAN.  It  has 
potential  for  growth  in  terms  of  versatility  and  problem 
size,  and  lends  itself  to  tailoring  to  suit  a  particular 
class  of  problem.  MDOT  could  be  coupled  with  a  microcomputer 
CAD  package  to  close  the  design  loop  in  the  computer. 

MDOT  has  been  validated  by  tests  on  a  number  of  problems, 
both  constrained  and  unconstrained.  Its  performance  is  good, 
and  could  be  made  better  through  refinement  of  the 
algorithms.  Specific  modifications  might  involve  the  one¬ 
dimensional  search  routines  and  the  frequency  of  gradient 
calculations . 

As  microcomputers  continue  to  become  more  commonplace  and 
their  capabilities  continue  to  improve,  emphasis  will  shift 
away  from  the  mainframes  for  the  solution  of  problems  which 
are  on  a  scale  easily  handled  by  smaller  machines.  Software 
such  as  MDOT  will  both  accompany  and  encourage  this  shift. 


APPENDIX  A 


MDOT  USER  MANUAL 

1.  INTRODUCTION 

To  avoid  repetition,  references  to  material  presented  in 
the  preceding  chapters  are  made  in  this  appendix.  A  useful 
follow-on  project  would  be  to  assemble  a  user  manual  for 
MDOT  independent  of  the  background  and  developmental 
material  in  the  body  of  the  thesis.  Details  of  computer 
operation  have  not  been  included  here,  as  it  is  assumed  that 
the  user  is  either  familiar  with  the  machine  or  has  access 
to  the  operating  manual. 

MDOT  is  currently  available  only  on  magnetic  tape 
cartridge  for  use  in  Hewlett-Packard  series  80  computers.  If 
it  is  to  be  implemented  in  an  HP-85A,  the  machine  must  be 
configured  with  four  enhancements:  a  16  Kbyte  memory  exten¬ 
sion  module,  a  matrix  ROM,  advanced  programming  ROM,  and  a 
printer  /  plotter  ROM.  No  peripheral  devices  are  required, 
nor  is  extensive  programming. 

2.  PROBLEM  FORMULATION 

Formulation  of  a  well-posed  problem,  as  discussed  in 
chapter  II,  is  fundamental  to  the  satisfactory  performance 
of  an  optimization  program.  First,  the  design  variables  must 
be  identified.  These  are  the  parameters  of  the  problem  which 


PfW'PjV .*y  *yjif.*^  .v-^7  v./.-jr;  rr* 


?-v  a.v,"-  v^-  '.’•  '>'av-  '.v'A". 


the  optimizer  will  be  permitted  to  change  in  its  search  for 
the  best  design.  The  objective  must  then  be  a  function  of 
these  variables,  and  the  minimum  of  this  function  is  what 
the  optimizer  will  seek.  Constraints  may  be  imposed  on  the 
design  in  two  ways: 

1.  Upper  and/or  lower  bounds  (side  constraints)  may  be 
specified  for  any  of  the  design  variables. 

2.  General  inequality  constraints  may  be  expressed  as 
functions  of  the  design  variables. 

Side  constraints  are  explicitly  assigned  when  the  initial 
design  estimate  is  entered  into  the  problem  subprogram. 
Inequality  constraints  must  be  formulated  as  quantities 
which  are  to  be  less  than  or  equal  to  zero.  Care  should  be 
exercised  to  avoid  redundant  or  otherwise  unnecessary  con¬ 
straints.  In  MDOT,  an  unconstrained  problem  has  neither  side 
constraints  nor  inequality  constraints.  MDOT  has  no  pro¬ 
vision  for  equality  constrained  problems. 

The  objective  function  can  be  any  characteristic  of  the 
design  expressible  mathematically  in  terms  of  the  design 
variables.  It  is  important  to  keep  in  mind  that  it  is  the 
minimum  of  this  function  which  is  sought.  If  the  problem  is 
formulated  around  an  objective  which  is  to  be  maximized, 
then  it  must  be  entered  in  such  a  way  that  MDOT  will  seek  to 
minimize  the  negative  of  this  objective. 
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An  important  consideration  during  problem  formulation  is 
the  range  of  orders  of  magnitude  of  the  design  variables. 
Optimizer  performance  is  best  in  a  design  space  in  which  the 
contours  of  constant  objective  function  are  concentric  hy¬ 
perspheres,  such  that  a  given  change  in  one  variable  has  the 
same  effect  on  the  objective  as  an  equal  change  in  any  other 
variable  [Ref.  11,  p- 173-  In  practice,  this  is  approximated 
by  scaling  the  design  variables  such  that  they  are  all  of 
the  same  order  of  magnitude,  or  nearly  so.  Some  optimizers 
do  this  automatically,  MDOT  does  not. 

Selection  of  the  initial  design  point  from  which  to 
start  MDOT  will  affect  its  performance  and  problem  run  time. 
Any  available  information  which  will  improve  the  initial 
approximation  should  be  used.  If  a  constrained  problem  is 
being  entered,  a  check  should  be  made  to  ensure  the  initial 
design  falls  within  the  side  constraints.  Although  MDOT  is 
equipped  to  handle  initially  infeasible  designs,  convergence 
will  likely  be  more  rapid  if  the  initial  design  is  free  of 
violated  inequality  constraints.  In  some  cases,  initial 
feasibility  may  be  a  difficult  thing  to  build  into  the 
problem  formulation.  MDOT  will  display  the  results  of  the 
first  design  evaluation  and  indicate  whether  or  not  it  is 
feasible.  At  that  time,  the  user  may  elect  to  proceed  with 
the  optimization  or  edit  the  initial  design  and  restart  the 


program. 


3.  PROBLEM  ENTRY 

MDOT  problem  entry  is  accomplished  by  editing  the  PROB 
subprogram.  With  this  module  loaded  into  memory  and  listed 
on  the  CRT,  modifications  are  made  on  the  program  lines 
noted  below. 

Line  10: 

The  file  name  of  the  subprogram  is  changed  to  any 
name  up  to  6  characters  in  length,  except  any  of 
those  already  assigned  to  MDOT  files. 

Line  100: 

Just  after  the  word  DATA,  two  integers  are  added, 
separated  by  a  comma.  The  first  is  the  number  of 
design  variables,  (N1),  the  second  is  the  number 
of  inequality  constraints,  (N2).  For  unconstrained 
problems  N2  is  always  zero. 

Lines  201-210: 

Initial  values  of  the  design  variables  are  en¬ 
tered,  beginning  with  XI  on  line  201  and  contin¬ 
uing,  one  variable  per  line.  If  the  problem  is 
constrained,  each  initial  value  is  followed  by  the 
lower  and  upper  bounds  assigned  to  the  correspond¬ 
ing  variable.  If  no  bound  is  to  be  specified,  th'i 
field  is  filled  by  an  "N". 


Lines  231-398: 


These  lines  are  available  for  defining  expressions 
to  be  used  in  the  design  evaluation.  These  state¬ 
ments  will  be  executed  prior  to  each  objective  or 
constraint  computation.  Whenever  the  design  vari¬ 
ables  are  used  in  this,  and  the  remaining  sections 
of  the  subprogram,  they  are  expressed  as  XI, 
X2, . . .XO. 

Lines  400-459: 

These  lines  are  available  for  defining  the  object¬ 
ive  function,  which  must  be  assigned  the  variable 
name  F. 

Lines  500-9000 

These  lines  are  available  for  defining  the  in¬ 
equality  constraint  functions,  which  must  be  sub¬ 
scripted  variables  named  G(i),  i=1,N2. 

Depending  on  the  complexity  of  the  problem,  the  user  may 
elect  to  use  any  BASIC  programming  structures  in  this  sub¬ 
program.  As  examples,  F0R-NEXT  loops,  FUNCTIONS,  SUB¬ 
ROUTINES,  and  even  other  SUBPROGRAMS  could  be  used  in  the 
problem  formulation.  With  about  7.3  Kbytes  of  memory  avail¬ 
able  for  constrained  problem  entry,  and  twice  that  for 
unconstrained  problems,  there  is  space  available  for  con¬ 
siderable  programming.  Variable  names  used  may  be  any  except 
those  included  in  the  PR0B  nomenclature,  as  defined  in 


Appendix  B.  If  subscripted  variables  are  to  be  used,  a  DIM 
(dimension)  statement  must  be  inserted  near  the  beginning  of 
the  subprogram.  When  problem  entry  is  complete,  the  subpro¬ 
gram  is  stored  under  the  new  file  name. 

Listings  of  the  subprograms  created  for  the  test  cases 
of  Chapter  IV  are  included  at  the  end  of  Appendix  B,  under 
the  names  "BANANA"  and  "BEAM".  Comparison  of  these  listings 
to  that  of  PROB,  which  is  the  unedited  version,  will  help  to 
illustrate  the  problem  entry  procedure  for  both  constrained 
and  unconstrained  cases. 

4.  PROGRAM  EXECUTION 

MDOT  is  started  either  by  execution  of  LOAD  "Autost"  and 
RUN  commands,  or  by  powering  up  the  computer  after  inserting 
the  tape  cartridge,  as  explained  in  section  III.B.4.  After 
the  "welcome"  graphic,  the  user  is  queried  as  to  the  problem 
file  name,  progress  display  options  and  output  format 
desired.  Optimization  then  proceeds. 

The  user  may  choose  to  monitor  the  optimization  closely, 
or  perhaps  not  at  all.  For  example,  if  a  constrained  problem 
is  being  run  for  the  first  time,  the  user  may  want  to  check 
for  rapidly  changing  design  variables  so  that  the  option  of 
editing  and  restarting  might  be  exercised.  On  the  other 
hand,  if  a  problem  known  to  be  well-behaved  is  being  re-run 
with  relatively  minor  changes,  the  user  may  elect  to 
"ignore"  MDOT  until  the  solution  is  obtained. 


Any  of  the  values  carried  by  the  DEFALT  subprogram,  as 
defined  in  the  listing  in  Appendix  B,  may  be  changed.  They 
may  be  permanently  modified  by  editing  and  re-storing 
DEFALT.  Alternatively,  a  default  value  may  be  changed  during 
program  execution,  whenever  the  edit  option  is  invoked,  by 
reassigning  its  value  from  the  keyboard.  If  the  program  is 
restarted,  though,  DEFALT  is  recalled. 


APPENDIX  B 


ANNOTATED  PROGRAM  LISTINGS 

In  an  effort  to  conserve  memory,  MDOT  was  coded  without 
remark  statements.  The  program  listings  are  therefore  pre¬ 
sented  in  an  annotated  format  to  aid  in  following  the  logic 
flowpaths.  Nomenclature  common  to  two  or  more  modules  is 
defined  first,  then  the  additional  nomenclature  unique  to 
each  module  immediately  precedes  the  applicable  segment  of 
the  code  listing.  Numbers  in  parentheses  in  the  function 
descriptions  refer  to  line  numbers  in  the  associated  program 
segment.  Parentheses  following  a  variable  name  in  the  nomen¬ 
clature  lists  indicate  vectors,  while  parentheses  enclosing 
a  comma  indicate  two-dimensional  arrays. 


COMMON  NOMENCLATURE 

A1()  Addresses  in  G()  of  violated  and  active  constraints 

A2()  Values  of  constraints  identified  in  A 1  ( ) 

C1()  Vector  of  integer  default  values 

C2()  Vector  of  non-integer  default  values 

D()  Current  objective  function  gradient 

D1()  Previous  iteration  objective  function  gradient 

F  Current  value  of  objective  function 

FO  Initial  value  of  objective  function 
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FI  Previous  iteration  value  of  objective  function 

G()  Vector  of  current  constraint  values 

G1(,)  Gradients  of  violated  and  active  constraints 

H(,)  Approximation  to  the  inverse  of  the  Hessian 

L()  Vector  of  design  variable  lower  bounds 

N1  Number  of  design  variables 

N 2  Number  of  inequality  constraints 

N3  Number  of  currently  violated  constraints 

N4  Number  of  currently  active  constraints 

P$  Problem  subprogram  file  name  character  string 

P0-P3  Polynomial  approximation  coefficients 

Q1  Iteration  counter 

Q2  Function  evaluation  counter 

Q3  Convergence  counter 

Q4  Convergence  counter 

R1  Golden  section  ratio 

R2  Golden  section  ratio 

SO  Search  direction  vector 

X()  Current  design 

X0()  Initial  design 

X 1 ( )  Previous  iteration  design 

X9()  Working  vector  of  perturbed  design  variables 
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PROGRAM  FLAGS 


Unless  otherwise  noted,  the  flags  listed  are  set  to  a 
value  of  1  when  the  associated  condition  exists. 

FLAG  CONDITION 

1  Maximum  number  of  iterations  has  been  exceeded 

2  Number  of  violated  constraints  greater  than  N1+2 

3  Search  vector  components  are  all  essentially  zero 

4  No  appreciable  move  parameter  can  be  found 

5  Convergence  to  optimum  has  occurred 

6  Variable  metric  algorithm  is  to  be  restarted 

7-8  Unassigned 

9  Subprogram  is  to  be  exited  without  execution 

10  Termination  has  occurred 

11  Progress  option  select:  1  =  data  display 

2  =  graphic  display 
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Module : 


1 


Autost 


Calls:  Autost  calls  all  other  modules  of  MDOT  into 

memory,  then  scratches  those  not  needed  for  the 
type  of  problem  (constrained  or  unconstrained) 

f 

being  run.  For  information  transfer,  Autost 

calls:  PROB  (edited),  UCONT,  CCONT,  LOGO,  and  s 

DEFALT.  | 


Function:  MDOT  main  calling  program  (10-370). 

Interactive  problem  initiation  (1000-1200). 
Program  loading  (2000-2070). 

Initial  design  display  generation  (3000-3130). 
Design/default  editing  (4000-4260). 


Nomenclature : 

A$  Interactive  query  response  (string) 

S$&A$  Concatenated  string  to  call  appropriate  optimizer 
D$  Display  string  ("FEASIBLE"  or  "INFEASIBLE") 

B$  Display  string  ("C1("  or  "C2(") 

NO  Working  variable  for  editing  default  values 

N  Interactive  query  response  (numeric) 

K0  If  =  0,  indicates  first  problem  of  run 


I 


10  !  Autost. . . . 

20  OPTION  BASE  1 

30  COM  , INTEGER  N1 , N2, N3, N4, G1 , Q2, Q3, 04, Cl < 10 >  , SHORT 

C2 ( 20 ) 

40  DIM  S'i [43,  A*  113,  D$ C  1 0  3 ,  B$ 1 3  3 
50  SHORT  XO  < 1 0 ) , FO , G  < 50 ),L(10),U(10), NO 
60  INTEGER  N,I,KO 
70  SM"CONT" 

30  GOSUB  2000 

90  CALL  "LOGO"  <2  SCRATCHSUB  "LOGO" 

100  CALL  "DEFALT" 

110  Q1,Q2,Q3,Q4,N3,N4,K0=O 
120  FOR  1=1  TO  11 
130  CFLAG  I 
140  NEXT  I 

150  WAIT  1000  <2  GOSUB  1000 
160  CLEAR  @  CALL  Pi  <  0  ) 

170  REDIM  X0(N1) 

ISO  IF  N2=0  THEN  210 

190  REDIM  G(N2) ,L(N1 ) ,U(N1 ) 

200  GOTO  220 

210  MAT  G=2ER(1)@  MAT  L=G@  MAT  U=L 
220  CALL  Pi  (  1,X0< ),FO,G< ),L( ),U( )  > 

230  GOSUB  3000 

240  IF  N2=0  THEN  Ai="U"  ELSE  Ai="C" 

250  CALL  Ai&Si  (  X0< > , FO, G< > , L< > , U< >  ) 

260  CFLAG  10  S  CFLAG  9 

270  CLEAR  0  DISP  S  DISP  "Select  option..,"  @  DISP  @  DISP  " 

1)  EDIT/RESTART" 

230  DISP  ©  DISP  "  2)  NEW  PROBLEM"  @  DISP  S  DISP  "  3)  E 

X I T" 

290  INPUT  N 

300  ON  N  GOTO  320,310,360 
310  SCRATCHSUB  Pi  ©  GOTO  100 
320  GOSUB  4000 

330  CALL  Pi  (  1,X0( ),FO,G( ),L( >,U< )  ) 

340  GOSUB  3000 
350  GOTO  250 
360  CLEAR 
370  END 


1000  IF  K0>0  THEN  1070 
1010  K0=1 

1020  CLEAR  0  DISP  0  DISP  "Select  problem  type  for  this  ru 
n  of  MDOT : "  0  DISP 

1030  DISP  "  1)  Constrained"  0  DISP  0  DISP  "  2)  Uncon 

strained"  0  INPUT  N 
1040  ON  N  GOTO  1050,1060 

1050  SCRATCHSUB  "UCONT"  0  SCRATCHSUB  "VMSRCH"  0  SCRATCHSUB  " 
NEWH"  0  GOTO  1070 

1060  SCRATCHSUB  "CCONT"  0  SCRATCHSUB  "FDSRCH"  0  SCRATCHSUB  " 
ACON"  0  SCRATCHSUB  "DIRECT" 

1070  CLEAR  0  DISP  0  DISP  "Have  you  created  ?<  stored  your  pr 
ob lem?" 

1030  DISP  0  DISP  "  (V  or  N;  ENTER)"  0  INPUT  At- 

1090  IF  A*="N"  THEN  1120 
1100  IF  A*="Y"  THEN  1160 
1110  GOTO  1000 

1120  CLEAR  0  DISP  0  DISP  "  Please  refer  to  the"  0  DISP 

1130  DISP  "  MDOT  USER  MANUAL"  0  DISP 

1140  DISP  "  for  instructions..." 

11 50  GOTO  370 

1160  DISP  0  DISP  "Enter  problem  subprogram  file  name. 

.."0  INPUT  P* 

1170  CLEAR  0  DISP  0  DISP  "Select  progress  option..."  0  DISP 
0  DISP  "  1)  GRAPHIC  DISPLAY" 

.  1130  DISP  0  DISP  "  2)  DATA  DISPLAY"  0  INPUT  N0  CFLAG  11 

1190  IF  N=2  THEN  SFLAG  11 
1200  RETURN 

2000  SFLAG  9  0  DISP  "Loading  MDOT..." 

2010  CALL  "DEFALT"  0  CALL  "CCONT"  0  CALL  "UCONT"  0  CALL  "ACO 
N" 

2020  DISP  0  DISP  "Still  loadinq  MDOT..." 

2030  CALL  "GRAD"  0  CALL  "DIRECT"  0  CALL  "FDSRCH"  0  CALL  "VMS 
RCH  " 

2040  DISP  0  DISP  "Almost  finished..." 

2050  CALL  "NEWH"  0  CALL  "CONV"  0  CALL  "PROG"  0  CALL  "TERM" 
2060  CFLAG  9  0  CLEAR 
2070  RETURN 
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3000  CLEAR  3  DISP  USING  3010  ;  Pt 

3010  IMAGE  3/ ,  13X , 6A 

3020  DISP  USING  3030  ;  N1,N2 

3030  IMAGE  1/,4X,2D,"  Design  Var iab les" , / , 4X , 2D, "  Inequality 
Constraints" 

3040  DISP  USING  3050  ;  F0 

3050  IMAGE  2/ , IX, "Initial  Designs  F  =  ",K 
3060  IF  N2=0  THEN  3090 

3070  IF  AMAX ( G) <C2( 3)  THEN  Dt=" FEASIBLE"  ELSE  Dt=" INFEASIBLE 

II 

30S0  DISP  USING  "  1 IX , 10A"  ;  D* 

3090  ON  KEY#  1,"EDIT"  GOSUB  4000 
3100  ON  KEY#  4, "CONTINUE"  GOTO  3130 
3110  KEY  LABEL 
3120  GOTO  3120 
3130  RETURN 

4000  CLEAR  a  DISP  a  DISP  "Select  editing  option..."  a  DISP  0 
DISP  "  1)  DESIGN  VARIABLES" 

4010  DISP  0  DISP  "  2)  DEFAULT  VALUES"  0  INPUT  N 

4020  ON  N  GOTO  4030,4140 

4030  CLEAR  0  DISP  "Enter  the  address  in  X  of  the  variable 

to  be  changed"  0  INPUT  I 

4040  DISP  USING  4050  ;  I,X0(I) 

4050  IMAGE  2/, IX, "current  value:  " , "X  < " , K, " )  =  ",K 
4060  DISP  0  DISP  "Enter  the  new  value..."  a  INPUT  X0(I) 

4070  CLEAR  0  DISP  0  DISP  "Editing  complete  ?"  a  DISP  " 

(Y  or  N;  ENTER)"  0  INPUT  A* 

4080  IF  A$="N"  THEN  4000 
4090  IF  A*="Y"  THEN  4110 
4100  GOTO  4070 
4110  Q1,Q2=0 

4120  CALL  PI  (  2, X0( ),F0,G( ),L( ),U( )  ) 

4130  GOTO  230 

4140  CLEAR  0  DISP  0  DISP  "Select  default  array..."  0  DISP  0 
DISP  "  1 )  CIO  integers" 

4150  DISP  0  DISP  "  2)  C2(>  reals"  0  INPUT  N 

4160  CLEAR  0  DISP  "Enter  the  address  in  C  of  the  variable 

to  be  changed"  0  INPUT  I 

4170  IF  N=1  THEN  N0=C1(I)  ELSE  N0=C2(I) 

41S0  IF  N= 1  THEN  B*="C1<"  ELSE  B*="C2<" 

4190  DISP  USING  4200  ;  B*,I,N0 

4200  IMAGE  2/, IX, "current  value:  ",3A,K," )  =  " , K 
4210  DISP  0  DISP  "Enter  the  new  value..."  0  INPUT  N0 
4220  IF  N=1  THEN  C1<I)=N0  ELSE  C2(I)=N0 

4230  CLEAR  0  DISP  0  DISP  "Editinq  complete  ?"  0  DIS®  "  ( 

Y  or  N;  ENTER)"  0  INPUT  A* 

4240  IF  A*="N"  THEN  4000 
4250  IF  A*="Y"  THEN  4260 
4260  RETURN 


Function:  Generate  MDOT  "welcome"  graphic  display  on  CRT. 
Scratched  from  memory  upon  execution. 


1  SUB  "LOGO" 

10  PEN  1  2  GCLEAR 
20  SCALE  0,10,0,10 
30  PEN  -1  a  GCLEAR  7 
40  PEN  1  a  GCLEAR  3 
50  CSI2E  30  a  LORG  5 
60  MOVE  5,5  a  LABEL 
70  CSI2E  6  a  MOVE  5,. 
30  WAIT  1000  a  PEN  1 
90  CSI2E  12  a  LORG  2 
100  MOVE  1,3  a  LABEL 
110  MOVE  1,6  a  LABEL 
120  MOVE  1,4  a  LABEL 
130  MOVE  1,23  LABEL 
140  CSI2E  3 
150  MOVE  2,3  a  LABEL 
160  MOVE  2,6  a  LABEL 
170  MOVE  2,4  a  LABEL 
130  MOVE  2,2  a  LABEL 
190  SUBEND 
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Module : 


DEFALT 


Called  by:  Autost 

Function:  Set  MDOT  working  parameters  to  their  default  val¬ 

ues.  DEFALT  is  called  each  time  a  new  problem  is 
executed . 


Nomenclature : 

C 1 ( 1 )  Maximum  number  of  iterations 

Cl (2)  Consecutive  iterations  for  convergence  criteria 
Cl (3)  Variable  metric  method  select:  0  =  DFP,  1  =  BFGS 
Cl (4)-C2( 10)  Unassigned 

C2(1)  Finite  difference  perturbation  factor 

C2(2)  Minimum  absolute  finite  difference  step 

C2(3)  Violated  constraint  criterion  (tolerance) 

C2(4)  Active  constraint  criterion  (thickness) 

C2(5)  Push-off  factor  multiplier  (theta  zero) 

C2(6)  Maximum  value  of  push-off  factor 

C2(7)  Factor  used  in  DIRECT  when  infeasible 

C2(8)  Factor  used  in  step  length  estimate  based  on  F 

C2(9)  Factor  used  in  step  length  estimate  based  on  X 

C2(10)  F  convergence  criterion  (relative) 

C2( 11)  F  convergence  criterion  (absolute) 

C2(12)  Defined  zero 

C2 ( 13)  Epsilon,  used  to  prevent  division  by  zero 

C2(14)-C2(20)  Unassigned 


y 

y 

y 


10  SUE  "DEFALT" 

20  OPTION  BASE  1 

30  COM  PIC  63  , INTEGER  N1 , N2, N3, N4, Q1 , Q2, Q3, Q4, Cl ( 10 )  , SHORT 
C2 ( 20 ) 

40  MAT  C1=2ER(10) 

42  MAT  C2=ZER(20) 

44  Cl ( 1 ) =20  !  max  #  iterations 
46  Cl(2)=2  !  consec.  conv. 

43  Cl(3)=l  w  in  H  update 

50  C2(l)=.01  !  fin.  diff.  mult. 

52  C2 < 2 ) = . 00 1  !  min.  fin.  step 
54  C2(3)=.004  !  const,  viol. 

56  C2 ( 4 ) =-. 1  !  active  const. 

53  C2(5)=l  !  push-off  mult. 

60  C2<6)=50  !  max.  push-off 

62  C2(7)=100000  !  Phi  (DIRECT) 

64  02(3)=. 1  !  obj.  mult. 

66  C2(9)=.l  !  des.  var.  mult. 

63  C2( 10)=.001  !  min  rel  F 
70  C2(ll)=.001  !  min  abs  F 
72  C2( 12) =.001  !  zero 

74  C2(13)=.0001  \  epsilon 
100  SUBEND 
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Module : 
Called  by: 


PROB 


User 


Function:  Provide  the  user  with  a  skeleton  subprogram  into 
which  the  optimization  problem  is  entered  by 
editing. 


Module:  PROB  (edited) 

Called  by:  Autost,  GRAD,  FDSRCH,  VMSRCH 

Function:  Input  the  number  of  design  variables  and  the 

number  of  inequality  constraints  (90-100). 

Input  the  intial  design  and  side  constraints 
(130-210) . 

Evaluate  the  objective  function  and  inequality 
constraints  (230-9000). 

Nomenclature: 

K1  Flag  to  indicate  first  call  or  subsequent  call 

L$  String  used  in  lower  bound  input 

U$  String  used  in  upper  bound  input 

X1-X0  Design  variable  names  used  in  problem  input 


>: 


i 

£ 


62 


10  SUB  "PROB"  ( K 1 , X  < ) , F , G ( ) , L ( ) , U ( ) ) 

20  OPTION  BASE  1 

30  COM  P*C6I)  ,  INTEGER  N1 ,  N2,  N3»  N4,  Q1 ,  02,  03,  04,  Cl  ( 10 ) 
C2 ( 20 ) 

35  SHORT  X1,X2,X3,X4,X5,X6,X7,X3,X9,X0 
40  DIM  L*C63,U*C63 
50  IF  K1 >1  THEN  220 
90  READ  N1,N2 
100  DATA 

105  IF  K1=0  THEN  SUBEXIT 
130  FOR  1=1  TO  N1 
140  READ  X ( I > 

142  IF  N2=0  THEN  170 
144  READ  L*»U* 

ISO  IF  L>=,,N"  THEN  L(I)=-1.E99  ELSE  L  (  I  >  =VAL  ( L*  ) 

160  IF  U*="N"  THEN  U(I)  =  1.E99  ELSE  U<  I  )=VALdJt  ) 

170  NEXT  I 

201  DATA 

202  DATA 

203  DATA 

204  DATA 

205  DATA 

206  DATA 

207  DATA 
203  DATA 

209  DATA 

210  DATA 

220  GOSUB  9010 

230  !  User-defined  expressions 
399  !  Objective  function 

499  02=02+1 

500  !  CONSTRAINTS 
9000  SUBEND 

9010  X1=X(1)  @  IF  N1 =1  THEN  RETURN 

9020  X2=X ( 2 )  ®  IF  Nl=2  THEN  RETURN 

9030  X3=X (3)  ®  IF  Nl=3  THEN  RETURN 

9040  X4=X ( 4 )  S  IF  Nl=4  THEN  RETURN 

9050  X5=X ( 5 )  @  IF  Nl=5  THEN  RETURN 

9060  X6=X ( 6 )  6  IF  N1 =6  THEN  RETURN 

9070  X7=X ( 7 )  ®  IF  Nl=7  THEN  RETURN 

9030  XS=X<3)  ®  IF  Nl=3  THEN  RETURN 

9090  X9=X ( 9 )  ®  IF  N1 =9  THEN  RETURN 

9100  X0=X(10)  @  RETURN 


Module : 


CCONT 


Called  by: 
Calls: 

Function : 


Autost 

TERM,  PROG,  ACON,  GRAD,  DIRECT,  FDSRCH,  CONV 

Control  constrained  optimization  by  the  method  of 
feasible  directions  (Fig.  2). 
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10  SUB  "CCONT"  (X0( ),F0,G( ),L(),U( )) 

20  OPTION  EASE  1 

30  COM  PIC 6 3  , INTEGER  N1 , N2, N3, N4, Q1 , Q2, Q3, Q4, Cl ( 10 )  , SHORT 
C2 ( 20 ) 

40  SHORT  X 1  ( 1 0 ) ,  X  (  1 0 ) ,  D ( 1 0 ) ,  G 1  < 12, 1 0 ) , A2 ( 1 2 ) , S ( 10) 

50  SHORT  F,B0,A0,F1 
60  INTEGER  Al(12) 

70  IF  FLAG ( 9 )  THEN  SUBEXIT 

30  REDIM  X1(N1),D(N1),G1(N1+2,N1),S(N1) 

90  A0,F1,B0=0  ©  F=F0 
100  MAT  S=ZER(N1>©  MAT  X=X0 
110  MAT  Al=ZER(Nl+2)©  MAT  A2=A1 
120  G 1 =0 1 + 1 

130  IF  Q 1 < C 1 ( 1 )  THEN  160 

140  SFLAG  1  ©  CALL  "TERM"  (  X(),F,G<)  )  ©  SUBEXIT 
150  SUBEXIT 

160  CALL  "PROG"  (  X0< ) , X  <  ) , F/F0, F  )  ©  IF  FLAG ( 9 )  THEN  SUBEXI 
T 

170  CALL  "ACQN"  (  G< ) , A1 < > , A2< )  )  ©  IF  FLAG ( 9 )  THEN  SUEEXIT 
180  IF  FLAG (2)  THEN  CALL  "TERM"  (  X(),F,G()  ) 

190  IF  FLAG( 10)  THEN  SUBEXIT 

200  IF  N3+N4#0  THEN  REDIM  G1(N3+N4,N1) 

210  CALL  "GRAD"  (  X ( ) , F , D <  > , G <  > , A 1 ( > , A2 ( ) , G 1 ( , )  )  ©  IF  FLAG ( 

9)  THEN  SUBEXIT 

220  IF  N3+N4#0  THEN  240 

230  MAT  S=-D©  GOTO  360 

240  CALL  "DIRECT"  (  D< > ,G1 < , > , A1 < > , A2< > ,S< > ,B0  >  @  IF  FLAG ( 9 
)  THEN  SUBEXIT 
250  CFLAG  3 

260  IF  MAXAB(S) >C2< 12)  AND  ABS(B0) >C2( 12)  THEN  360 
270  SFLAG  3 

280  IF  N3#0  THEN  CALL  "TERM"  (  X(),F,G()  ) 

290  IF  FLAG (10)  THEN  SUBEXIT 

300  IF  MAXAB(  A2X=C2(  12)  THEN  CALL  "TERM"  (  X(),F,G()  ) 

310  IF  FLAG( 10)  THEN  SUBEXIT 

320  IF  C2(4) >=-C2( 12)  THEN  CALL  "TERM"  <  X(),F,G()  ) 

330  IF  FLAG(10)  OR  FLAG ( 9 )  THEN  SUBEXIT 
340  C2 ( 4 )=C2 ( 4 ) /3 
350  GOTO  170 

360  MAT  S=( 1/MAXAB(S) )*S©  CFLAG  4 
370  MAT  X1=X©  F 1 =F 

380  CALL  "FDSRCH"  (  X ( ) , L ( ) , U ( ) , F , D ( ) , G ( ) , A 1 ( ) , A2 ( ) , G 1 ( , ) , S ( 
),A0  ) 

390  IF  FLAG(9)  THEN  SUBEXIT 
400  IF  FLAG ( 4 )  THEN  280 
410  CFLAG  5 

420  MAT  X=< 1)*X1+(A0)*S 
430  CALL  Ft  (  2,XO,F,G(>  ) 

440  CALL  "CONV"  (  X ( ) , X 1 ( ) , F, FI  )  ©  IF  FLAG ( 9 )  THEN  SUBEXIT 
450  IF  FLAG ( 5 )  THEN  CALL  "TERM"  (  X(),F,G()  ) 

460  IF  FLAG(10)  OR  FLAG ( 9 )  THEN  SUBEXIT  ELSE  GOTO  120 
470  SUBEND 


SS 


Module: 


UCONT 


Called  by:  Autost 

Calls:  GRAD,  PROG,  TERM,  VMSRCH,  CONV,  NEWH 

Function:  Control  unconstrained  optimization  by  the 
variable  metric  method  (Fig.  1). 

Nomenclature : 

K  Indicates  an  iteration  in  which  no  move  parameter 

to  improve  the  design  was  found 
A  One-dimensional  search  move  parameter 


1  0  SUE  ” UCONT "  (  XU  <  )  ,  F0 ,  G  <  >  ,  L  (  )  ,  U  (  )  ) 

20  OPTION  BASE  1 

30  COM  Pi C  6 3  ,  I NTEGER  N 1 , N2 , N3 , N4 , Q 1 , Q2 , Q3 , Q4, C 1(10)  , SHORT 
C2 <  20 ) 

40  SHORT  X< 10) , XI < 10) ,D( 10) ,01 < 10) ,S( 10) , H( 10, 10) 

50  SHORT  F, FI , A 
60  INTEGER  K 

70  IF  FLAG<9)  THEN  SUBEXIT 
S0  REDIM  X(N1) , X1(N1) , D ( N 1 ) , D 1 ( N 1 ) 

90  MAT  X=X0@  F=F0  ©  K=0 

100  CALL  "GRAD"  (  X ( ) , F, D( ) , G( )  ) 

110  MAT  H=IDN(  N1 , N1 ) 

120  Q 1 =Q 1+1 

130  CALL  "PROG"  (  X  0 ( ) , X ( ) , F / F 0 , F  )  ©  IF  FLAG  <  9 )  THEN  SUEEXI 
T 

140  IF  QKCl(l)  THEN  170 

150  SFLAG  1  ©  CALL  "TERM"  <  X(),F,G()  ) 

160  IF  FLAG (10)  OR  FLAG <  9 )  THEN  SUBEXIT 
170  MAT  S=H*D 
130  MAT  S=-S 

190  MAT  S=( 1/MAXAB(S) >  *S 
200  MAT  X1=X©  MAT  01=0©  F1=F 

210  CALL  " VMSRCH"  (  X ( ) , F , 0 ( ) , S ( ) , A  )  ©  IF  FLAG ( 9 )  THEN  SUEE 
X  IT 

220  IF  A >C2<12)  THEN  270 
230  K=K+1 

240  IF  K<2  THEN  110 

250  SFLAG  4  ©  CALL  "TERM"  (  X(),F,G()  ) 

260  IF  FLAG<10)  OR  FLAG(9)  THEN  SUBEXIT 
270  MAT  X=< 1 )*X+(A)*S 

230  CALL  "CQNV"  (  X < ) , XI ( ) , F, Fl  )  ©  IF  FLAG ( 9 )  THEN  SUBEXIT 
290  IF  FLAG  ( 5 )  THEN  CALL  "TERM"  (  XO,F,G()  ) 

300  IF  FLAG( 10)  OR  FLAG<  9)  THEN  SUBEXIT 
310  IF  FLAG( 6)  THEN  110 

320  CALL  "GRAD"  <  X ( ) , F , 0 <  ) , G ( )  )  ©  IF  FLAG < 9 )  THEN  SUBEXIT 
330  CALL  "NEWH"  (  X ( ) , X 1 ( ) , D ( ) , 01 ( ) , H ( , )  )  ©  IF  FLAG  <  9 )  THEN 
SUBEXIT 
340  GOTO  120 
350  SUBEND 
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Module:  ACON 


Called  by:  CCONT 


Function:  Determine  the  number  of  currently  violated  (N3), 
and  active  (N4),  inequality  constraints. 

Construct  the  A1()  vector  of  the  addresses  in  G() 
of  the  active/violated  constraint  set. 

Construct  the  A2()  vector  of  the  current  values 
of  this  set. 

Violated  constraint  information  is  stored  in  the 
first  N3  rows  of  A1()  and  A2(),  active  constraint 
information  in  the  last  N4  rows. 


V'.V.O.N* 


10  SUB  "ACON"  <G< ) , A1 ( ) , A2< ) ) 

20  OPTION  BASE  1 

30  COM  PIC 63  ,1 NTEGER  N 1 , N2 , N3 , N4 , Q 1 , 02 , Q3 , Q4 , 
C2 <  2® ) 

40  ON  KEY#  1  GOTO  300 

50  IF  FLAG ( 9 )  THEN  SUBEXIT 

60  CFLAG  2  2  N3,N4=0 

70  FOR  1=1  TO  N2 

30  IF  G ( I ) >=C2 ( 3 )  THEN  110 

90  IF  G ( I ) >=C2 ( 4 )  THEN  N4=N4+1 

100  GOTO  120 

110  N3=N3+1 

120  NEXT  I 

130  IF  N3<  =N 1 +2  THEN  150 

140  SFLAG  2  2  SUBEXIT 

150  IF  N3+N4=0  THEN  SUBEXIT 

160  MAT  A1=2ER ( N3+N4 ) 2  MAT  A2=A1 

170  03,04=1 

130  FOR  1=1  TO  N2 

190  IF  G < I ) >=C2 ( 3 )  THEN  250 

200  IF  G < I ) < C2 < 4 )  THEN  280 

210  A1 ( N3+04 ) = I 

220  A2 ( N3+04 ) =G ( I ) 

230  04=04+1 
240  GOTO  280 
250  A 1 ( 03 ) = I 
260  A2 <  03 ) =G < I ) 

270  03=03+1 

280  NEXT  I 

290  GOTO  310 

300  SFLAG  9  2  SUBEXIT 

310  SUBEND 


Module : 


GRAD 


Called  by:  CCONT,  UCONT 
Calls:  PROB  (edited) 

Function:  Calculate  the  gradient  of  the  objective  function 
by  first  forward  finite  difference  approximation. 
In  constrained  optimization,  calculate  gradients 
of  the  active/violated  inequality  constraint  set. 

Nomenclature: 

F2 


N 


Intermediate  function  evaluation 
Design  variable  perturbation 


1 0  SUB  "GRAB"  (X< ),F,D(),G( ),A1< ),A2( ),G1(,  )  ) 

2®  OPTION  BASE  1 

3®  COM  P*Z£>1  , INTEGER  N1 , N2« N3, N4, Q1 , Q2, Q3,Q4,C1 < 10) 
C2 ( 20 ) 

40  SHORT  F2,X?<10),N 

50  ON  KEY#  1  GOTO  210 

60  IF  FLAG < 9 )  THEN  SUBEXIT 

70  FOR  1=1  TO  N1 

S0  MAT  X9=X 

90  N=C2 ( 1 )*ABS<  X9( I ) ) 

100  IF  N<C2 ( 2 )  THEN  N=C2<2) 

110  X9<I )=X9<I )+N 

120  CALL  P*  <  2,X9< ),F2,G< )  ) 

130  IF  N2=0  THEN  IS® 

14®  IF  N3+N4=0  THEN  IS® 

15®  FOR  J=1  TO  N3+N4 

16®  G1(.J,  I  )  =  < G ( A 1  < •  J )  )-A2(  J)  )/N 

17®  NEXT  J 

13®  D< I )=<F2-F)/N 

19®  NEXT  I 

200  GOTO  22® 

210  SFLAG  9  @  SUBEXIT 
22®  SUBEND 


Module : 


DIRECT 


Called  by:  CCONT 

Function:  Solve  the  direction-finding  subproblem  in  the 
method  of  feasible  directions. 

Calculate  constraint  push-off  factors  (170-200). 
Initialize  working  arrays  for  currently  feasible 
(2000-2080),  or  infeasible  (1000-1090)  designs. 
Determine  the  direction  vector  SO,  and  the  para¬ 
meter  B0  (300-700). 


Nomenclature : 

A( , )  Working  array  constructed  from  G1(,)  and  TO 
B( , )  Working  array  initialized  as  -ATA 

B0  Kuhn-Tucker  parameter 

B2  Intermediate  element  value  used  in  pivoting 

B3  Intermediate  element  value  used  in  pivoting 

B9  Intermediate  variable  used  in  pivoting 

C()  Working  vector 

D0()  Working  vector  initialized  as  D() 

G0(,)  Working  array  initialized  as  G1(,) 

G9()  Working  vector  used  in  constructing  G0(,) 

I9()  Working  vector  of  element  indices 

J9  Working  integer  scaler 

K9  Working  integer  scaler 

N9  Dimension  of  active/violated  constraint  set 


Working  vector 

Vector  of  constraint  push-off  factors 


10  SUB  "DIRECT"  (D(),G1(, >,A1<>,A2(>,S< ),B0> 

20  OPTION  BASE  1 

30  COM  P*C6I  , INTEGER  N1 , N2, N3, N4, Q1 , Q2, Q3, Q4 , Cl < 10 )  , SHORT 
C2 ( 20 ) 

40  SHORT  B2, B3,B9, A< 13, 11 ) ,T(12) *B< 13» 13) ,P( 11 ) ,U< 13) ,C( 13) , 
YU1),D0<  10),  00(12,  10),  09(10) 

50  INTEGER  J9 , K9 , I 9 ( 1 3 ) , N9 
60  ON  KEY#  1  GOTO  650 
70  IF  FLAG (9)  THEN  SUBEXIT 
80  N9=N3+N4 

90  REDIM  T(N9),P(N1+1),Y(N1+1) 

100  FOR  1=1  TO  N9 

110  T(I)=C2(5)*<1-A2(I)/C2(4) )~2 

120  IF  T ( I ) >C2 ( 6 )  THEN  T ( I ) =C2 ( 6 ) 

130  NEXT  I 

140  MAT  P=ZER(Nl+l)g  MAT  B0=D 
150  MAT  D0=( 1/MAXAB(D0) )*D0 
160  MAT  G0=G 1 
170  FOR  1=1  TO  N9 
130  MAT  G9=G0 (  I ,  ) 

190  MAT  G9= ( 1 /MAX AB ( G9 ) ) *G9 
200  MAT  G0(I,)=G9 
210  NEXT  I 

220  IF  N3>0  THEN  GOSUB  1000  ELSE  GOSUB  2000 
230  MAT  B=A*TRN ( A ) 

240  MAT  B=-B 
250  N=UBND(B,1) 

260  MAT  I9=ZER ( N ) 

270  A9, K9=0 

280  FOR  1=1  TO  N 

290  IF  C(I)>=0  THEN  330 

300  B9=C(I)/B(I, I) 

310  IF  B9<=A9  THEN  330 
320  K9=I  S  A9=B9 
330  NEXT  I 

340  IF  K9=0  THEN  540 
350  J9=I9(K9 ) 

360  I9(K9)=0 

370  IF  J9=0  THEN  I9(K9)=K9 
380  B2*B<K9,K9) 

390  FOR  1=1  TO  N 
400  B(K9,I)=B(K9, I)/B2 
410  NEXT  I 


420  C(K9)=A9 
430  B  <  K9 , K9 ) = 1 / B2 
440  FOR  1=1  TO  N 
450  IF  I=K9  THEN  520 
460  B3=B<I,K9) 

470  B  < I , K9 ) =0 

480  FOR  J=1  TO  N 

490  B< I, J)=B( I, J)-B3*B<K9, J) 

500  NEXT  J 

510  C<I)=C<I)-B3*A9 

520  NEXT  I 

530  GOTO  270 

540  FOR  1=1  TO  N 

550  U<I)=0 

560  J9=I9< I ) 

570  IF  J9<  =0  THEN  590 
530  U  < I ) =C  <  J9 ) 

590  NEXT  I 
600  MAT  Y=TRN<A)*U 
610  MAT  Y=P— Y 
620  MAT  S=Y< IsNl ) 

630  B0=Y <  Nl  +  1 ) 

640  GOTO  660 

650  3FLAG  9  @  SUBEXIT 

660  SUBEND 

1000  MAT  C=ZER<N9)@  MAT  U=ZER<N9) 

1010  MAT  A=CQN ( N9 , N 1 + 1 ) 

1020  MAT  A< , 1:N1 )=G0 
1030  MAT  A( ,  Nl  +  1 ) =T 
1040  MAT  P<1:N1)=D0 
1050  MAT  P=-P 
1060  P  ( N 1  + 1 ) =C2  <  7 ) 

1070  MAT  C=A*P 
1080  MAT  C=-C 
1090  RETURN 

2000  MAT  C=ZER<N9+1 )2  MAT  U=ZER(N9+1) 
2010  MAT  A*CGN<N9+1 ,  Nl  +  1 ) 

2020  MAT  A  < 1 : N9 , 1 s  N 1 ) =00 
2030  MAT  A<1:N9,N1+1)=T 
2040  MAT  A  <  N9+ 1 , 1 : N 1 ) =D0 
2050  P<N1+1 ) =1 
2060  MAT  C=A<*N1+1) 

2070  MAT  C=-C 
2080  RETURN 


Module:  FDSRCH 

Called  by:  CCONT 
Calls:  PROB  (edited) 

Function:  Perform  one-dimensional  search  for  constrained 

optimization  in  the  method  of  feasible 
directions . 

Estimate  initial  search  move  parameter  (1000- 
1180) . 

Check  for  side  constraint  violations  (2000-2050). 
Establish  bounds  on  solution,  feasible  (3000- 
3340)  or  infeasible  (4000-4360). 

Refine  solution  by  polynomial  approximation, 
feasible  (5000-5240)  or  infeasible  (6000-6580). 

Nomenclature : 

A  Move  parameter 

A0()  Working  vector  of  move  parameters 

A()  Working  vector  of  move  parameters 

A1  Initial  A  based  on  change  in  objective  function 

A2  Initial  A  based  on  attaining  feasibility 

B(,)  Array  of  constraint  values  during  search 

DO  Dot  product  of  D()  and  S()  or  G 1 ( i , )  and  S() 

GOO  Working  vector  of  constraint  gradients 
M()  Working  vector  of  maximum  constraint  values 

Y()  Working  vector  of  objective  function  values 
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10  SUB  "FDSRCH"  ( X ( ) , L ( ) , U( ) , F, D( ) , G( ) , A1 < ) , A2( ) , G1 < , ) , S< ) , A 
) 

20  OPTION  BASE  1 

30  COM  Ft Z 63  ,1 NTEOER  N1 ,  N2,  N3,  N4 ,  Q 1 ,  Q2, Q3 ,04,01(1®)  , SHORT 
C2 ( 2® ) 

40  SHORT  X9< 10 ) , R1 , R2, A1 , A2, A ( 4 ) , Y < 4 ) , M< 4 ) , B ( 4, 50 ) , D®, 0© ( 10 ) 
,P®,P1,P2,B®,A®(51) 

5©  INTEGER  I,K 

60  ON  KEY#  1  GOTO  21© 

7©  IF  FLAG  (?)  THEN  SUBEXIT 

80  REBIM  X9(N1),B(4, N2 ) , G® ( N 1 ) , A® ( N2+ 1 ) 

9©  R1=(3-SQR(5) )/2  @  R2=2-R1 
10©  A , A 1 , A2 , A3 , A ( 1 ) =®  ©  CFLAG  4 

110  MAT  M=ZER<4>©  MAT  X9=X©  Y(1)=F  ©  MAT  B<1,)=0©  M(1)=AMAX( 
G) 

12©  GOSUB  1000 

13©  A ( 2 ) =A  ©  GOSUB  2®0© 

140  IF  N3=©  THEN  GOSUB  3®®®  ELSE  GOSUB  4000 
15©  IF  N3=0  THEN  GOSUB  500©  ELSE  GOSUB  6©0© 

160  A=A(K) 

17©  IF  A >C2 (12)  THEN  19© 

13©  SFLAG  4  ©  SUBEXIT 
19©  F=Y<K)  @  MAT  G®B(K, ) 

20©  GOTO  22© 

21©  SFLAG  9  0  SUBEXIT 
22©  SUBEND 
I©©©  GOSUB  105© 

1010  IF  N3#0  THEN  1©3© 

102©  A=A1  @  RETURN 

103©  IF  A2>2*A1  THEN  A=2»A1  ELSE  A=MAX(A1,A2) 

104©  RETURN 
105©  D©=DOT ( D, S ) 

106©  A1=C2(8)*ABS(F) /ABS(D® ) 

107©  FOR  1=1  TO  N1 

108©  A=C2 ( 9 ) * ABS ( X ( I ) ) /ABS(S( I ) ) 

109©  IF  A< A1  THEN  A1=A 
110©  NEXT  I 

111©  IF  N3=0  THEN  RETURN 
112©  FOR  1=1  TO  N3 
113©  MAT  G®=G1 ( I , ) 

114©  D©=DOT(G0,S) 

115©  A=-(A2( I )/D0) 

1160  IF  A>A2  THEN  A2=A 
1170  NEXT  I 
1130  RETURN 

2000  MAT  X=< 1 )*X9+(A>*S 

201©  FOR  1=1  TO  N1 

2020  IF  X(IXL(I)  THEN  X(I)=L(I) 

2030  IF  X(I)>IJ(I)  THEN  X(I)=U(I) 

204©  NEXT  I 
205©  RETURN 
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300<j*  M<  1  >=AMAX(G) 

3010  CALL  Pi  <  2,X(),F,G(>  ) 
3020  M  <  2 )=AMAX ( G ) 

3030  Y  <  2 ) =F 

3040  MAT  E( 2, )=G 

3050  IF  Y < 2 >  >Y ( 1 >  THEN  3070 

3060  IF  M ( 2 ) <  =C2 (12)  THEN  3150 

3070  A(3)=A< 1 )+Rl*< A(2)-A< 1 ) ) 

3080  A=A ( 3 ) 

3090  GOSUB  2000 

3100  CALL  Pi  <  2,X(),FtG()  ) 

3110  M(3)=AMAX(G) 

3120  Y(3)=F 
3130  MAT  B(3» )=G 
3140  RETURN 
3150  A (  3 ) =A  <  2 ) 

3160  Y(3)=Y(2) 

3170  M< 3)=M( 2) 

3130  MAT  G=B ( 2, ) 

3190  MAT  B( 3, )=G 

3200  A  <  2 )  =  ( 1 +R2 ) * A ( 3 ) -R2* A ( 1 ) 

3210  A=A ( 2 ) 

3220  GOSUB  2000 

3230  CALL  Pi  (  2,X(),F,G<>  ) 

3240  Y ( 2 ) =F 

3250  M ( 2 ) =AMAX ( G ) 

3260  MAT  B<2t)»G 
3270  IF  Y  ( 2 )  >Y  ( 3 )  THEN  RETURN 
3280  IF  M ( 2 ) >C2 (12)  THEN  RETURN 
3290  A  < 1 ) =A  <  3 ) 

3300  Y ( 1 ) =Y ( 3 ) 

3310  M ( 1 )=M(3) 

3320  MAT  G=B<3,) 

3330  MAT  B<1, )=G 
3340  GOTO  3150 


j  .i .  —w .  mm .  — jwi'mwiw  ■. 


T 


■'.*> vV y TO w>".  ’-'WA ^ .'•  ^'AV- 


4000  M( 1 )=AMAX  <G) 

4010  CALL  PI  (  2,X(),F,G()  ) 

4020  Y(2)=F 

4030  M( 2) =AMAX ( G ) 

4040  MAT  B ( 2, )=G 
4050  IF  11(2)  >11(1)  THEN  4030 
4060  IF  M ( 2 ) >0  THEN  4160 
4070  IF  Y ( 2 ) < Y ( 1 )  THEN  4160 
4080  A(3)=A( 1 )+Rl*(A(2)-A( 1 ) > 
409 0  A=A ( 3 ) 

4100  GOSUB  2000 

4110  CALL  P*  <  2,X(),F,G()  ) 

4120  Y ( 3 ) =F 

4130  M(3)=AMAX(G) 

4140  MAT  B(3»)*G 
4150  RETURN 
4160  A<3)=A<2) 

4170  Y ( 3 ) =Y ( 2) 

4180  MAT  G=B ( 2, ) 

4190  MAT  B ( 3, )=G 
4200  M ( 3 ) =M( 2 ) 

4210  A<2)=< 1+R2)*A<3)-R2*A< 1 ) 
4220  A=A ( 2 ) 

4230  GOSUB  2000 

4240  CALL  P*  (  2,X(),F,G()  ) 

4250  Y(2)=F 

4260  M ( 2 ) =AMAX ( G ) 

4270  MAT  B(2,)=G 
4280  IF  M( 2) >M( 3)  THEN  RETURN 
4290  IF  M ( 2 ) >0  THEN  4310 
4300  IF  Y < 2 ) >Y ( 3 )  THEN  RETURN 
4310  A( 1 )=A( 3) 

4320  Y ( 1 ) = Y ( 3 ) 

4330  M( 1 )=M(3) 

4340  MAT  G=B ( 3, ) 

4350  MAT  B( 1, )=G 
4360  GOTO  4160 
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5000  P2=  (  ( Y  ( 3 )  -  Y  (  1  )  )  /  (  A  (  3  )  -  A  <  1  )  )  -  <  Y  (  2 )  -  Y  (  1  )  )  /  <  A  <  2 )  -  A  <  1  )  )  >  /  <  A 
(3)-A(2) ) 

5010  P 1 = ( Y ( 2 ) — Y (1))/(A(2) -A ( 1 ) )-P2*<A< 1 )+A<2) ) 

5020  A0 ( 1 ) =- ( P 1 / ( 2*P2 ) ) 

5030  FOR  1=1  TO  N2 

5040  P2= ( <B<3, I )-B( 1, I ) )/(A(3)-A( 1 ) )-<B<2, I >-B< 1, I ) )/(A(2)-A 
( 1 ) ) ) / ( A ( 3 ) -A ( 2 ) ) 

5050  P 1  =  ( B  ( 2 , I )-B( 1, I ) )/(A(2)-A(l ) )-P2*<A< 1 )+A<2) ) 

5060  P0=B < 1 , I ) -P 1 *A  < 1 > -P2*A  < 1 ) A2 
5070  B0=P1'"2-4*P0*P2 
5030  IF  B0>0  THEN  5110 
5090  A0( 1+1 )=-(P0/Pl ) 

5100  GOTO  5130 

5110  A0 ( I + 1 ) =MAX ( -P 1 +SGR ( B0 ) , -P 1 -SGR ( B0 ) ) 

5120  A0< 1+1 )=A0< 1+1 )/<2*P2> 

5130  NEXT  I 

5140  A(4)=AMIN(A0) 

5150  A=A  <  4 ) 

5160  GOSUB  2000 

5170  CALL  P*  (  2,X(),F,G()  ) 

5180  Y <  4 ) =F  @  MAT  B<4, >=G 
5190  M ( 4 ) =AMAX ( G ) 

5200  K= 1 

5210  FOR  1=2  TO  4 

5220  IF  YdXY(K)  AND  M(IX=C2(12)  THEN  K=I 
5230  NEXT  I 
5240  RETURN 


.'*v  1  >  t  7  ^7  7; 


y  TrTTTPFT?  -■'  -OW1-  / 


6000  MAT  A0=ZER  <  N3+ 1  ) 

6005  IF  AM  IN  ( M ) >0  THEN  6360 

6010  P2=< ( Y ( 3 ) -Y  < 1 ) ) / ( A ( 3 ) -A ( 1 ) ) - ( Y ( 2 ) - Y  < 1 ) ) / <  A  <  2 ) -A  < 1 ) ) ) / ( A 
(3)-A(2) ) 

6020  P 1 = ( Y ( 2 ) — Y ( 1 ) )/( A(2)-A< 1 ) )-P2*< A< 1 )+A<2) ) 

6030  A0(1)=-(P1/(2*P2) ) 

6040  FOR  1=1  TO  N3 
6050  GOSUB  6490 
6060  IF  B0>0  THEN  6090 
6070  A0 ( 1+ 1 ) =- ( P0/P1 ) 

6080  GOTO  6110 

6090  A0 ( I + 1 ) =M I N ( -P 1 +SQR ( B0 ) , -P 1 -SQR ( B0 ) ) 

6100  A0 ( I + 1 ) =A0 ( I + 1 ) / ( 2*P2 ) 

6110  NEXT  I 

6120  A(4)=AMAX( A0) 

6130  MAT  A0=ZER ( N2-N3 ) 0  Jl=l 

6140  FOR  1=1  TO  N2 

6150  FOR  J=1  TO  N3 

6160  IF  I=A1 ( J )  THEN  6240 

6170  NEXT  J 

6130  GOSUB  6540 

6190  IF  B0>0  THEN  6220 

6200  A0(.J1  )=-(P0/Pl  ) 

6210  Jl-Jl+1  ®  GOTO  6240 

6220  A0  ( J 1 )  =M  I N  (  -P 1 +SQR  ( B0 )  ,  -P 1  -SQR  ( B0 )  ) 

6230  A0( J1 )=A0( J1 >/<2*P2) 

6235  J1=J1+1 
6240  NEXT  I 
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•  .  *  •  '  *  .•  ./  .  *  A  *  ,  *  *«  W  * 


6250  A<4)=MIN(A<4) , AMIN(A0> ) 

6260  A=A ( 4 ) 

6270  GOSUB  2000 

6280  CALL  P*  (  2,XOtFtG<>  ) 

6290  Y ( 4 ) =F  ®  MAT  B<4,)=G 
6300  M(4)=AMAX(G) 

6310  K=1 

6320  FOR  1=2  TO  4 

6330  IF  Y <  I ) < Y < K )  AND  M(IX=C2(4)  THEN  K=I 
6340  NEXT  I 
6350  RETURN 

6360  P2=<  <  M  ( 3 )  —M  ( 1 )  )  /  <  A  ( 3 )  — A  ( 1 ) )  —  ( M  ( 2 )  — M  <  1 )  )  /  <  A  <  2 )  -  A  (  1  )  )  )  /  ( A 
( 3 ) -A ( 2 )  ) 

6370  P1=(M(2)-M( 1 ) )/ (A(2)-A( 1 ) >-P2*(A( 1 )+A(2) ) 

6380  A ( 4 ) =- ( P 1 / <  2*P2 ) ) 

6390  IF  A ( 3 ) < A < 4 )  AND  A(4XA(2)  THEN  6420 
6400  M ( 4 ) =2*M ( 2 ) 

6410  GOTO  6460 
6420  A=A ( 4 ) 

6430  GOSUB  2000 

6440  CALL  P*  (  2 , X ( ) , F , G  < )  ) 

6450  M(4)=AMAX(G)  @  Y<4)=F  @  MAT  B(4,)=G 
6460  M=AMIN<M) 

6470  K=AM INROW 
6480  RETURN 

6490  P2=<  <B<3, A1 ( I ) ) -B< 1 , A1 ( I ) ) ) / ( A(3 >-A ( 1 ) ) -( B ( 2, A1 ( I ) ) -B ( 1 
,  A 1  ( I )  )  )  /  ( A  <  2 )  -  A  (  1  )  )  )  /  ( A  ( 3 )  -A  <  2 )  ) 

6500  P1=(B<2, A1 < I ) )-B( 1, A1 (I ) ) )/(A(2)-A( 1 ) )-P2*<A( 1 )+A(2) ) 
6510  P0=B<  1 ,  A1  ( I )  )-Pl*A(  1  )-P2*A(  1  )  --2 
6520  B0=P1A2-4*P0*P2 
6530  RETURN 

6540  P2=( (B(3, I )-B< 1 , I ) )/ ( A(3)-A( 1 ) )-(B(2, I )-B( 1 , I ) )/ (A (2) -A 
( 1 ) ) )/<A(3)-A<2) ) 

6550  P1=(B(2, I )-B< 1 , I ) )/<A<2)-A( 1 ) )-P2*( A( 1 )+A(2) ) 

6560  P0=B  ( 1 ,  I )  -P 1  *A  ( 1 )  -P2*A  <  1 )  ''2 
6570  B0=P1''2-4*P0*P2 
6580  RETURN 


Module : 


VMSRCH 


Called  by:  UCONT 


Calls : 


PROB  (edited) 


Function:  Perform  unconstrained  one-dimensional  search  in 

variable  metric  method. 

Establish  bounds  on  the  solution  (100-350). 

Refine  the  solution  by  polynomial  approximation 
(360-520). 


Nomenclature : 

A0  Move  parameter 

A( )  Working  vector  of  move  parameters 

A1()  Working  vector  of  move  parameters 

DO  Dot  product  of  D()  and  SO 

F9  Working  value  of  objective  function 

Y()  Working  vector  of  objective  function  values 

Y1()  Working  vector  of  objective  function  values 


rtr* 


10  SUB  " VMSRCH"  (X( ),F,D( ),S( ),A0) 

20  OPTION  BASE  1 

30  COM  P*  C  6 1  , I NTEGER  N 1 , N2 , N3 , N4 , Q 1 , Q2 , Q3 , Q4 , C 1 ( 1 © )  , SHORT 
C2 ( 20 ) 

40  SHORT  DQ, A1 < 1 1 ) , X9( IQ) , Y (4 > , A< 4 ) , R1 , R2, PI , P2,P3, Y1 ( 4 ) , F9 , 

B 

50  ON  KEY#  1  GOTO  540 

60  IF  FLAG < 9)  THEN  SUEEXIT 

70  R 1  =  ( 3-SGR ( 5 ) ) / 2  ®  R2= < 1+SQR< 5) > /2 

80  MAT  A=ZER (4)0  MAT  Y=A 

90  Y ( 1 ) =F  S  F9=0 

100  D0=DOT (D,S) 

110  MAT  A1=ZER ( Nl+1 ) 

120  MAT  X9=X 

130  A1 ( 1 )=C2(3)*ABS(F)/ABS(D0> 

140  FOR  1=1  TO  N1 

150  A1 ( 1  +  1 )=C2<9)*ABS( X9( I ) ) / ABS ( S < I ) ) 

160  NEXT  I 

170  A(2)=AMIN(A1 ) 

180  MAT  X9=< 1 )*X9+<A<2> ) *S 

190  CALL  P*  <  2, X9< > ,F9  >  ; 

200  Y ( 2 ) =F9 

210  IF  Y ( 2 ) >Y ( 1 )  THEN  310 
220  A(3)=A(2)  ®  Y(3)=Y(2) 

230  A ( 2 )  =  ( 1 +R2 ) * A ( 3 ) -R2* A  < 1 > 

240  MAT  X9=X 

250  MAT  X9=< 1)*X9+<A<2> )*S  j 

260  CALL  Pt-  (  2,X9<>,F9  )  ! 

270  Y ( 2 ) =F9  .  ; 

280  IF  Y <  2 ) >Y ( 3 )  THEN  360 
290  A  < 1 ) =A ( 3 )  @  Y ( 1 ) =Y  <  3 ) 

300  GOTO  220 

310  A ( 3 ) =R1*A ( 2 )  j 

320  MAT  X9=X  J 

330  MAT  X9=( 1 )*X9+< A(3) >*S  j 

V 

t 

1 


! 


> 

f 

I 
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340  CALL  P*  (  2tX9()tF9  ) 

350  V  <  3 ) =F9 

360  P3= ( Y ( 3 ) -V ( 1 ) >*<A(2>-A< 1 ) )/<A(3)-A( 1 >-<Y<2)-Y< 1 >  >*< A  (3)- 
A< 1 )  )  ) 

370  P3=P3/ ( A < 2 ) -A ( 1 )  +  ( A<3)-A<2> )*D0> / ( (A(2)-A( 1 ) )*< A(3)-A( 1 ) 
)*<A<3)-A<2> ) ) 

330  P2=< ( Y  ( 2 )  —  Y  ( 1 ) )  /  <  A  ( 2 )  — A  ( 1 ) > -D0 > / < A< 2 ) -A < 1 ) ) -P3* < 2*A ( 1 )+A 
(2)  ) 

390  P 1 =D0-2*P2*A ( 1 ) -3*P3*A < 1) A2 

400  B=P2'S2-3*P1*P3 

410  IF  B>=0  THEN  430 

420  J=3  Q  GOTO  490 

430  A ( 4 ) = ( -P2+SGR ( B ) ) / < 3*P3 ) 

440  MAT  X9=X 

450  MAT  X9=< 1 )*X9+<A<4> )  *S 
460  CALL  P*  (  2, X9( ) , F9  ) 

470  Y ( 4 ) =F9 
430  J=4 

490  MAT  Y1=Y ( 1 s  J) 

500  F=AMIN(Y1) 

510  K= AM INROW 
520  A0=A(K) 

530  GOTO  550 

540  SFLAG  9  @  SUBEXIT 

550  SUBEND 


Module : 


NEWH 


Called  by:  UCONT 

Function:  Update  the  approximation  to  the  inverse  of  the 
Hessian  matrix  used  in  determining  the  search 
direction  in  the  variable  metric  method. 


Nomenclature : 


D0(,) 

Update  matrix 

PO 

Working  vector 

initialized 

33  XO-XIO 

S 

Dot  product  of 

PO  and  YO 

TO 

Working  vector 

TIC) 

Working  vector 

T2(  ,  ) 

Working  array 

T3(,) 

Working  array 

Y() 

Working  vector 

initialized 

as  DO-DIO 
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1 0  SUE  "NEWH"  < X<  )  , XI < ) ,D( ) ,D1 ( ) ,H( ,  )  ) 

20  OPTION  EASE  1 

30  COM  PTC 61  ,  I NTEGER  N 1 , N2, N3 , N4 , Q 1 , Q2 , Q3 , Q4 ,C1(10)  , SHORT 
C2 ( 20 ) 

40  SHORT  D0( 10, 10) , Y( 10) ,P( 10) ,T1 ( 10) ,T2( 10, 10) ,T( 1 ) ,S,T3< 10 
,  10) 

50  ON  KEY#  1  GOTO  230 
60  IF  FLAG < 9 )  THEN  SUEEXIT 

70  REDIM  D0(N1,N1),Y(N1),P(N1),T1(N1),T2(N1,N1),T3(N1,N1) 

80  MAT  P=X-X1 
90  MAT  Y=D-D1 
100  S=DQT ( P, Y) 

110  MAT  T1=H*Y 

120  MAT  T2=T1*TRN(T1 ) 

130  MAT  T=TRN(Y)*T1 
140  MAT  D0=P*TRN(P) 

150  MAT  D0=< (S+T< 1 )*C1 (3) )/Sr2)*B0 
160  MAT  D0=( 1 )*D0+< (Cl (3)-l )/T( 1 ) )*T2 
170  MAT  T2=T1*TRN(P) 

180  MAT  T3=P*TRN(T1) 

190  MAT  T2=T2+T3 

200  MAT  D0= ( 1 ) *D0+  <-(Cl(3)/S) ) *T2 

210  MAT  H=H+D0 

220  GOTO  240 

230  SFLAG  9  ©  SUBEXIT 

240  SUBEND 


Module 


CONV 


Called  by:  CCONT,  UCONT 

Function:  Determine  whether  the  design  has  converged  to  the 

optimum  in  the  last  iteration. 

Update  convergence  criteria  based  on  iteration 
history. 
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iu 


Jriv  u. 


10  SUB  "CONV"  <X< ),X1< ) ,F,F1) 

20  OPTION  BASE  1 

30  COM  ,  INTEGER  N1 ,  N2,  N3,  N4,  Q1 ,  02,  03,  04,  Cl  <  10 )  , SHORT 

C2  <  20 ) 

40  SHORT  X9< 10) 

50  IF  FLAG < 9 )  THEN  SUBEXIT 
60  CFLAG  6 

70  C2 < 8 )  =  < C2 < 8 ) +ABS <  < F 1 -F ) /FI ) )/2 
80  MAT  X9=X1-X@  MAT  X9=X9/X1 
90  C2 <9 )  =  ( C2 (9 ) +MAXAB ( X9 ) ) /2 

100  IF  ABS<F1-FKMIN<C2<13),C2<11)*ABS<F1>)  THEN  03=03+1  ELS 
E  03=0 

110  IF  ABS(F1-F)/MAX(ABS(F1) , .00001)<C2( 10)  THEN  04=04+1  ELS 
E  Q4=0 

120  IF  MAX  < Q3 , 04 ) >=C 1 < 2 )  THEN  SFLAG  5 
130  IF  MAX (Q3, Q4 ) >0  THEN  SFLAG  6 
140  SUBEND 


10  SUB  "PROG"  <XO<),XO,Y,F) 

20  OPTION  BASE  1 

30  COM  PH 6 3  , INTEGER  N1 ,  N2,  N3,  N4,  Q1 ,  Q2,  Q3,  Q4,  Cl  ( 10 )  ,  SHORT 
C2  <  20 ) 

40  ON  KEY#  1 , " INTRPT"  GOTO  300 

50  IF  FLAG<9)  THEN  SUBEXIT 

60  IF  FLAG (11)  THEN  230 

70  IF  Q1 >1  THEN  ISO 

SO  GCLEAR  @  PEN  1  @  LORG  5 

90  SCALE  -2, 21 ,-3. 2, 3. 2 

100  XAXIS  0, 1,0,20 

110  YAXIS  0, 1,-3, 3 

120  FOR  I=-3  TO  3 

130  MOVE  -.5,1  0  LABEL  I 

140  NEXT  I 

150  MOVE  10,3  ®  LABEL  "Iteration  History" 

160  MOVE  10,2.7  a  LABEL  PI 

170  MOVE  -3.1,5  @  LABEL  "K1  to  interrupt" 

130  FOR  1=1  TO  N1 

190  MOVE  Q1 , X  < I ) /XO < I )  a  LABEL  I 

200  NEXT  I 

210  MOVE  Q1,Y  0  LABEL  "f" 

220  SUBEXIT 

230  IF  01  >1  THEN  270 

240  CLEAR  a  DISP  USING  "13X,6A"  ;  P$ 

250  DISP  "Iteration  Objective  Function" 

260  KEY  LABEL 

270  DISP  USING  2S0  ;  Q1,F 

280  IMAGE  3X , 2D, 13X , K 

290  GOTO  310 

300  SFLAG  9  a  SUBEXIT 

310  SUBEND 
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Module:  TERM 

Called  by:  CCONT,  UCONT 

Function:  Generate  output  of  optimization  results. 
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10  SUB  "TERM"  <X<  > ,F,G<  >  > 

20  OPTION  BASE  1 

30  COM  P*C63  , INTEGER  N1 , N2, N3, N4, G1 , 02, 03, 04, Cl ( 10 )  , SHORT 
C2 ( 20 ) 

40  IF  FLAG ( 9 )  THEN  SUBEXIT 
50  GOSUB  1000 

60  IF  FLAG ( 1 )  THEN  GOSUB  2000 

70  IF  FLAG (2)  THEN  GOSUB  3000 

80  IF  FLAG (3)  THEN  GOSUB  4000 

90  IF  FLAG<4)  THEN  GOSUB  5000 

100  IF  FLAG(5)  THEN  GOSUB  6000 

110  IF  FLAG (11)  THEN  150 

120  PRINT  3  PRINT  3  PRINT 

130  GRAPH  3  COPY  3  SUBEXIT 

140  PRINT  3  PRINT  3  COPY 

150  SUBEND 

1000  PRINT  3  PRINT 

1010  PRINT  USING  "  13X ,  6A"  ;  PI- 

1020  PRINT  USING  1030  ;  F 

1030  IMAGE  1 / , 3X , "OPTIMUM  =  " , K 

1040  FOR  1=1  TO  N1 

1050  PRINT  USING  1060  ;  I,X(I) 

1060  IMAGE  1/, 10X, "X( ",K, " )  =  ",K 
1070  NEXT  I 

1 030  IF  N2=0  THEN  1130 

1090  FOR  1=1  TO  N2 

1100  PRINT  USING  1110  ;  I , G ( I ) 

1110  IMAGE  1/, 10X, "G< ",K, " )  =  " , K 
1120  NEXT  I 

1130  PRINT  3  PRINT  3  PRINT  "Termination  based  on:"  3  PRINT 
1140  SFLAG  10 
1150  RETURN 

2000  PRINT  "Exceeded  max.  no.  of  iterations."  3  RETURN 
3000  PRINT  "Excessive  number  of  violated  constraints."  3 
RETURN 

4000  PRINT  "Failure  to  find  a  direction  to  improve  the  desi 
gn."  3  RETURN 

5000  PRINT  "Failure  to  find  a  move  parameter  to  improve  the 

design."  3  RETURN 

6000  PRINT  "Convergence"  3  RETURN 
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Listings  of  the  edited  versions  of  PROB  which  were  used 
for  the  test  cases  presented  in  Chapter  IV  are  included  as 
examples  of  MDOT  problem  input.  The  subprogram  used  to  enter 
the  unconstrained  test  problem  was  renamed  "BANANA",  while 
that  edited  for  the  constrained  problem  was  renamed  "BEAM". 


10  SUB  "BANANA"  ( ki  ,  X  (  ) ,  F , G ( ) , L (  )  ,  U (  )  ) 

20  OPTION  BASE  1 

30  COM  p* C 6 2  , INTEGER  N 1 , N2 , N3 , N4 , Q 1,02, Q3 , Q4 , C 1 < 1 0 )  , SHOP T 
C2<20) 

35  SHORT  XI,  X2,  X3, X4, X5, X6, X7, X8, X9, X0 
40  DIM  LiC6],U-^C63 
50  IF  K1>1  THEN  220 
90  READ  N1,N2 
100  DATA  2,0 

105  IF  K1=0  THEN  SUBEXIT 
130  FOR  1=1  TO  N1 
140  READ  X ( I ) 

142  IF  N2=0  THEN  170 
144  READ  LT, UT 

150  IF  LT="N"  THEN  L<I>  — 1.E99  ELSE  L ( I ) =VAL ( L* ) 

160  IF  UT="N"  THEN  U(I)  =  1.E99  ELSE  U  < I ) =VAL  < U T ) 

170  NEXT  I 

201  DATA  -1 

202  DATA  1.5 

203  DATA 

204  DATA 

205  DATA  • 

206  DATA 

207  DATA 

208  DATA 

209  DATA 

210  DATA 

220  GOSUE  9010 

230  !  User-def ined  expressions 

399  !  Objective  function 

400  F=10*Xl''4-20*Xl''-2*X2+10*X2"'"2+Xl"'2-2*Xl+5 

499  Q2=Q2+1 

500  !  CONSTRAINTS 
9000  SUBEND 

9010  X1=X(1)  @  IF  Nl=l  THEN  RETURN 

9020  X2=X (2)  @  IF  Nl=2  THEN  RETURN 

9030  X3=X<3)  S  IF  Nl=3  THEN  RETURN 

9040  X4=X(4)  @  IF  Nl=4  THEN  RETURN 

9050  X5=X (5)  @  IF  Nl=5  THEN  RETURN 

9060  X6=X ( 6 )  ©  IF  Nl=6  THEN  RETURN 

9070  X7=X ( 7 )  ©  IF  Nl=7  THEN  RETURN 

9080  X8=X ( 8 )  ©  IF  Nl=3  THEN  RETURN 

9090  X9=X ( 9)  ©  IF  Nl=9  THEN  RETURN 

9100  X0=X(10)  ©  RETURN 


* 
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10  SUB  "BEAM"  <K1,X<  >,F,G<  ),LO,U<  )> 

20  OPTION  BASE  1 

30  COM  P*C63  , INTEGER  N1 , N2, N3, N4, Q1 , 02, 03, 04, Cl < 10 ) 
C2 ( 20 ) 

33  SHORT  X1,X2,X3,X4,X5,X6,X7,X3,X9,X0 
40  DIM  L*C63,U*C63 
50  IF  K1 >1  THEN  220 
90  READ  N1,N2 
100  DATA  2,3 

105  IF  K1=0  THEN  SUBEXIT 

110  REDIM  L<N1),U(N1),X<N1),G(N2) 

130  FOR  1=1  TO  N1 
140  READ  X(I) 

142  IF  N2=0  THEN  170 
144  READ  L*,U* 

150  IF  L*="N"  THEN  L(I)=-1.E99  ELSE  L(I)=VAL(L*> 

160  IF  U*="N"  THEN  U(I)=1.E99  ELSE  U ( I ) =VAL ( U* ) 

170  NEXT  I 

201  DATA  3.5, .5,5 

202  DATA  16,1,20 

203  DATA 

204  DATA 

205  DATA 

206  DATA 

207  DATA 

208  DATA 

209  DATA 

210  DATA 

220  GOSUB  9010 

230  !  User-defined  expressions 

301  B=X 1 

302  H=X2 

399  !  Objective  function 

400  F=200*B*H 

499  02=02+ 1 

500  !  CONSTRAINTS 

501  G<1)=600/(B*H~2)-1 

502  G< 2) =10666. 7/ (B*H^3)-1 

503  G(3)=H/10-B 
9000  SUBEND 

9010  X1=X(1)  3  IF  Nl=l  THEN  RETURN 

9020  X2=X(2)  a  IF  Nl=2  THEN  RETURN 

9030  X3=X (3)  a  IF  Nl=3  THEN  RETURN 

9040  X4=X ( 4 )  8  IF  Nl=4  THEN  RETURN 

9050  X5=X <5 )  3  IF  Nl=5  THEN  RETURN 

9060  X6=X (6)  a  IF  Nl=6  THEN  RETURN 

9070  X7=X ( 7 )  a  IF  Nl=7  THEN  RETURN 

9080  X8=X ( 8 )  8  IF  Nl=8  THEN  RETURN 

9090  X9=X ( 9 )  a  IF  Nl=9  THEN  RETURN 

9100  X0=X(10)  a  RETURN 
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